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Abstract 

According to a theorem of Poincare, the solutions to differential equations are analytic 
functions of (and therefore have Taylor expansions in) the initial conditions and various 
parameters providing the right sides of the differential equations are analytic in the vari- 
ables, the time, and the parameters. We describe how these Taylor expansions may be 
obtained, to any desired order, by integration of what we call the complete variational 
equations. As illustrated in a Duffing equation stroboscopic map example, these Taylor ex- 
pansions, truncated at an appropriate order thereby providing polynomial approximations, 
can well reproduce the behavior (including infinite period doubling cascades and strange 
attractors) of the solutions of the underlying differential equations. 

1 Introduction 

In his prize-winning essay[l,2] and subsequent monumental work on celestial mechanics[3], 
Poincare established and exploited the fact that solutions to differential equations are 
frequently analytic functions of the initial conditions and various parameters (should any 
occur). This result is often referred to as Poincare analyticity or Poincare 's theorem on 
analyticity. 

"Corresponding author. 
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Specifically, consider any set of m first-order differential equations of the form 



Za = faizi,--- ,Zm;t;Xi,--- ,Xn), a = l,---,m. (1.1) 

Here t is the independent variable, the Za are dependent variables, and the A5 are possible 
parameters. Let the quantities z^ be initial conditions specified at some initial time t = t^, 

Za{t') = Zl (1.2) 

Then, under mild conditions imposed on the functions /„ that appear on the right side of 
(1.1) and thereby define the set of differential equations, there exists a unique solution 

Za{t) = ga{zl,- ■ ■ ,zl^;t^,t;\i,- ■ ■ ,Xn), a = l,m (1.3) 

of (1.1) with the property 

^a(i°) =5a(^?,--- ,^^;i°,i°;Ai,--- ,An) = ^°, a = l,m. (1.4) 

Now assume that the functions fa are analytic (within some domain) in the quantities 
Za, the time t, and the parameters A;,. Then, according to Poincare's Theorem, the solution 
given by (1.3) will be analytic (again within some domain) in the initial conditions z^, the 
times and t, and the parameters A;,. 

Poincare established this result on a case- by-case basis as needed using Cauchy 's method 
of majorants. It is now more commonly established in general using Picard iteration, and 
appears as background material in many standard texts on ordinary differential equations 
[4]. 

If the solution Za{t) is analytic in the initial conditions and the parameters Af,, then 
it is possible to expand it in the form of a Taylor series, with time-dependent coefficients, in 
the variables z^ and A^. The aim of this paper is to describe how these Taylor coefficients 
can be found as solutions to what we will call the complete variational equations. 

To aid further discussion, it is useful to also rephrase our goal in the context of maps. 
Suppose we rewrite the set of first-order differential equations (1.1) in the more compact 
vector form 

z = f{z;t;X). (1.5) 
Then, again using vector notation, their solution can be written in the form 

z(t)=5(zO;tO,i;A). (1.6) 

That is, the quantities z{t) at any time t are uniquely specified by the initial quantities z^ 
given at the initial time t^. 

We capitalize on this fact by introducing a slightly different notation. First, use f 
instead of t^ to denote the initial time. Similarly use to denote initial conditions by 
writing 

z' = z^ = z{f). (1.7) 
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Next, let be some final time, and define final conditions by writing 

zf = z{t^). (1.8) 

Then, with this notation, (1.6) can be rewritten in the form 

zf = g{z';t\tf;X). (1.9) 

We now view (1.9) as a map that sends the initial conditions to the final conditions 
z^. This map will be called the transfer map between the times f and , and will be 
denoted by the symbol Ai . What we have emphasized is that a set of first-order differential 
equations of the form (1.5) can be integrated to produce a transfer map Ai. We express 
the fact that ^A sends z^ to z-l" in symbols by writing 

zf = Mz\ (1.10) 

and illustrate this relation by the picture shown in Figure 1. We also note that M. is 
always invertible: Given z^ , t-^, and t*, we can always integrate (1.5) backward in time 
from the moment t = to the moment t = and thereby find the initial conditions z*. 
In the context of maps, our goal is to find a Taylor representation for M. If parameters 
are present, we may wish to have an expansion in some or all of them as well. 




Figure 1: The transfer map M. sends the initial conditions to the final conditions z-l" . 

The organization of this paper is as follows: Section 2 derives the complete variational 
equations without and with dependence on some or all parameters. Sections 3 and 4 
describe their solution using forward and backward integration. As an example. Section 5 
treats the Duffing equation and describes the properties of an associated stroboscopic map 
A4 . Section 6 sets up the complete variational equations for the Duffing equation, including 
some parameter dependence, and studies some of the properties of the map obtained by 
solving these variational equations numerically. There we will witness the remarkable 
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fact that a truncated Taylor map approximation to A4 can reproduce the infinite period- 
doubling Feigenbaum cascade and associated strange attractor exhibited by the exact M. 
Section 7 describes how the variational equations can be solved numerically. A final section 
provides a concluding summary. 

2 Complete Variational Equations 

This section derives the complete variational equations, first without parameter depen- 
dence, and then with parameter dependence. 

2.1 Case of No or Ignored Parameter Dependence 

Suppose the equations (1.1) do not depend on any parameters A;, or we do not wish to 
make expansions in them. We may then suppress their appearance to rewrite (1.1) in the 
form 

Za = fa{z,t), a = l,m. (2.1) 

Suppose that z'^{t) is some given design solution to these equations, and we wish to study 
solutions in the vicinity of this solution. That is, we wish to make expansions about this 
solution. Introduce deviation variables Ca by writing 

Za = zi + Ca. (2.2) 

Then the equations of motion (2.1) take the form 

zi + Ca = fa{z'' + C,t). (2.3) 

In accord with our hypothesis of analyticity, assume that the right side of (2.3) is analytic 
about z"^. Then we may write the relation 

Uz" + C, i) = /a(z^ t) + gaiz", t, C) (2.4) 
where each ga has a Taylor expansion of the form 

ga{z^t,0 = Y,mGr{0. (2.5) 

r 

Here the Gr{C) are the various monomials in the m variables Cfe labeled by an index r using 
some convenient labeling scheme, and the 5^ are (generally) time-dependent coefficients 
which we call forcing terms} By construction, all the monomials occurring in the right 
side of (2.5) have degree one or greater. We note that the Qait) are known once z'^{t) is 
given. 

^Here and in what follows the quantities ga are not to be confused with those appearing in (1.3). 
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By assumption, z'^ is a solution of (2.3) and therefore satisfies the relations 

zi = fa{z^t). (2.6) 
It follows that the deviation variables satisfy the equations of motion 

L=9a{z\t,Q = Y,gl{t)Gr{0. (2.7) 

r 

These equations are evidently generalizations of the usual first-degree (linear) variational 
equations, and will be called the complete variational equations. 

Consider the solution to the complete variational equations with initial conditions C,l 
specified at some initial time i*. Following Poincare, we expect that this solution will be 
an analytic function of the initial conditions C^- Also, since the right side of (2.7) vanishes 
when all Cb = [all the monomials Gr in (2.7) have degree one or greater], C,{t) = is a 
solution to (2.7). It follows that the solution to the complete variational equations has a 
Taylor expansion of the form 

ut) = Y.hmGr{C) (2.8) 

r 

where the h\{t) are functions to be determined, and again all the monomials Gr that occur 
have degree one or greater. When the quantities }fa{t) are evaluated at some final time t^ , 
(2.8) provides a representation of the transfer map M. about the design orbit in the Taylor 
form 

Ci = Ut^) = Y.^a{tf)Gr{C)- (2.9) 



2.2 Complete Veiriational Equations with Pcirameter Dependence 

What can be done if we desire to have an expansion in parameters as well? Suppose that 

there arc n such parameters, or that we wish to have expansions in n of them. The work of 
the previous section can be extended to handle this case by means of a simple trick: View 
the n parameters as additional variables, and "augment" the set of differential equations 
by additional differential equations that ensure these additional variables remain constant. 

In detail, suppose we label the parameters so that those in which we wish to have an 
expansion are Ai • • • A„. Introduce n additional variables Zm+i, ■ ■ ■ Z£ where i = m + n hy 
making the replacements 

\b Zm+b, b = l,n. (2-10) 
Next augment the equations (1.1) by n more of the form 

Za = 0, a = m + l,£. (2.11) 
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By this device we can rewrite the equations (1.1) in the form 

Za = fa{z,t), a = l,e (2.12) 

with the understanding that 

/„ = /a(z;i;y"^^), a = l,m, (2.13) 

where A"^^™ denotes the other remaining parameters, if any, and 

/a = 0, a = m + l,t (2.14) 

For the first m equations we impose, as before, the initial conditions 

Za{f) = zl a = l,m. (2.15) 

For the remaining equations we impose the initial conditions 

Za{f) = Xa-m, a = m + l,t (2.16) 

Note that the relations (2.14) then ensure the Za for a > m retain these values for aU t. 

To continue, let z'^{t) be some design solution. Then, by construction, we have the 
result 

zi{t) = Xtm = >^a-m, a = m + l,t (2.17) 
Again introduce deviation variables by writing 

Za = zi + Ca a = l,i. (2.18) 

Then the quantities Co for a > m will describe deviations in the parameter values. More- 
over, because we have assumed analyticity in the parameters as well, relations of the forms 
(2.4) and (2.5) will continue to hold except that the Gr{C,) are now the various monomials 
in the i variables C,},. Relations of the forms (2.6) and (2.7) will also hold with the provisos 
(2.13) and (2.14) and 

5^(i)=0, a = m + \,t (2.19) 

Therefore, we will only need to integrate the equations of the forms (2.6) and (2.7) for 
a < m. Finally, relations of the form (2.9) will continue to hold for a < m supplemented 
by the relations 

Cf = C, a = m + l,L (2.20) 

Since the Gr{C) now involve ^ variables, the relations of the form (2.9) will provide an 
expansion of the final quantities C^i (for a < m) in terms of the initial quantities (for 
a < m) and also the parameter deviations (^^ with a = m + l,i. 
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3 Solution of Complete Variational Equations Using For- 
ward Integration 

This section and the next describe two methods for the solution of the complete variational 
equations. This section describes the method that employs integration forward in time, 
and is the conceptually simpler of the two methods. 

3.1 Method of Forwcird Integration 

To determine the functions let us insert the expansion (2.8) into both sides of (2.7). 
With r" as a dummy index, the left side becomes the relation 

ia = Y.h'l'it)Gr»{C). (3.1) 

For the right side we find the intermediate result 

Y.gl{t)Gr{Q = Y.m GlY,hi{t)GAC),---Y.<^t)GAC^ . (3.2) 

r r \ r' r' / 

However, since the Gr are monomials, there are relations of the form 

Gr(Y. hi{t)GAC). • • • E hi{t)GAC)^ = E u;"{K)Gr'^{C), (3.3) 

and therefore the right side of (2.7) can be rewritten in the form 

Y^mGrio = j2J29im^iK)GAC). (3.4) 

r r" r 

The notation UJ!" (h^) is employed to indicate that these quantities might (at this stage of 
the argument) depend on all the with n ranging from 1 to m, and s ranging over all 
possible values. 

Now, in accord with (2.7), equate the right sides of (3.1) and (3.4) to obtain the relation 

j^/if (i)G,„(c) = Y.Yg:m;"{hi)GAC). (3.5) 

Since the monomials Gr"(C*) linearly independent, we must have the result 

h::'{t) = j29:{tM"{K). (3.6) 
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We have found a set of differential equations that must be satisfied by the h^^. Moreover, 
from (2.8) there is the relation 

Un=T.f'a{tlGr{C) = C. (3.7) 
r 

Thus, all the functions h'a(t) have a known value at the initial time f, and indeed are 
mostly initially zero. When the equations (3.6) are integrated forward from t = t'^iot = t^ 
to obtain the quantities /^^(t-^), the result is the transfer map M. about the design orbit in 
the Taylor form (2.9). 

Let us now examine the structure of this set of differential equations. A key observation 
is that the functions U^. (hf^) are universal. That is, as (3.3) indicates, they describe certain 
combinatorial properties of monomials. They depend only on the dimension m of the system 
under study, and are the same for all such systems. As (2.7) shows, what are peculiar to 
any given system are the forcing terms Qait)- 

3.2 Application of ForwEird Integration to the Two-veiriable Case 

To see what is going on in more detail, it is instructive to work out the first nontrivial case, 
that with m = 2. For two variables, all monomials in (Ci)C2) are of the form (Ci)"'^ (C2)''^- 
Here, to simplify notation, we have dropped the superscript i. Table 1 below shows a 
convenient way of labeling such monomials, and for this labeling we write 

G.(C) = (Ci)^nC2)^'^ (3.8) 

with 

h=h{r) and j2=32{r). (3.9) 



Table 1: A labeling scheme for monomials in two variables. 



r 


Ji 


J2 


D 


1 


1 





1 


2 





1 


1 


3 


2 





2 


4 


1 


1 


2 


5 





2 


2 


6 


3 





3 


7 


2 


1 


3 


8 


1 


2 


3 


9 





3 


3 
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Thus, for example, 

Gi = Ci, (3.10) 

G2 = C2, (3.11) 

G3 = Ci, (3.12) 

G4 = C1C2, (3.13) 

G5 = Cl etc. (3.14) 

For more detail about monomial labeling schemes, see Section 7. 

Let us now compute the first few {hfi)- Prom (3.3) and (3.10) we find the relation 

Gl(Ylh'lGr'{0,Y1^2Gr'{C?j = Ylh'lGr'{0 = J^Ul" Gr'' {()■ (3-15) 
It follows that there is the result 

C/["=/ir. (3.16) 
Similarly, from (3.3) and (3.11), we find the result 

U'2" = h'2- (3.17) 
Prom (3.3) and (3.12) we find the relation 

Gs fe^^'^-'(0,E^2G.'(c)') = fe^^'^'-'(o) 

\ r' r' / \ r' ) 

= hlh\GsiOGtiO = J2 ufGr'^iQ. (3.18) 

s,t r" 

Use of (3.18) and inspection of (3.10) through (3.14) yields the results 

= 0, (3.19) 

Ul = 0, (3.20) 

Ul = ih\f, (3.21) 

ul = 2h\hj, (3.22) 

Ul = {hlf. (3.23) 
Prom (3.3) and (3.13) we find the relation 



GAY,^i(GAC).Y.^iGAC) \ = [Y.k(gac)] (E^2 

\ r' r' / \ r' J \ r' 



Gr'iO 



J2 KhlGsiOGtiO = J2 UfGr'^iC). (3.24) 



9 



It follows that there are the results 



Ul = 0, 

ui = 0, 
u! = h\hi 



(3.25) 
(3.26) 
(3.27) 

(3.28) 
(3.29) 



= h\hl + hlh\, 
ul = hlhl. 

Finally, from (3.3) and (3.14), we find the results 



Ul = 0, 

Ul = 0, 
ul = {h\f, 
ul = 2hlhl 

Ul = [hlf. 



(3.30) 



(3.31) 

(3.32) 
(3.33) 
(3.34) 



Two features now become apparent. As in Table 1, let D{r) be the degree of the 
monomial with label r. Then, from the examples worked out, and quite generally from 
(3.3), we see that there is the relation 



It follows that the sum on the right side of (3.6) always terminates. Second, for the 
arguments h"^ possibly appearing in Ul (/?-fJ, wc see that there is the relation 



Taylor expansion (2.9) through terms of some degree D, it is only necessary to integrate a 
finite number of equations, and the right sides of these equations involve only the coefficients 
for this degree and lower. 

For example, to continue our discussion of the case of two variables, the equations (3.6) 
take the explicit form 



[/; = when D{r) > D{r"). 



(3.35) 




2 




(3.37) 



2 



h\{t) = Y.9mul = 9l{t)h\{t) + gl{t)h\{t) 



(3.38) 
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2 

hlit) = ^9[m? = 9l{t)hlit) + 9f{t)hlit), (3.39) 

r=l 
2 

h1{t) = Y.gimr = 9l{t)hlit) + gl{t)hlit), (3.40) 

r=l 
r=l 

= glimlit) + gf{t)hl{t) + gfrnhlit)? 

+gf{t)hl{t)hlit) + gmhl{t)]'', (3.41) 

hlit) = j^gim^r 

r=l 

= gl(t)h\{t) + 9lit)hl{t) + gl(t)[h\{t)f 

+52(i)/i}W/i2W+5i(t)[/i2(i)]', (3.42) 

h\{t) = Y.9imr 

r=l 

= glimiit) + gl{t)hi{t) + 2gf{t)h\{t)hl{t) 

+gt{t)[hl{t)hl{t) + hlmlit)] + 2gl{t)hl{t)hlit), (3.43) 

r=l 

+gl{t)[h\{t)hl{t) + hl{t)hl{t)] + 2glit)hl{t)hl{t), (3.44) 

hlit) = Y.9im'r 

r=l 

= gl{t)hUt) + gf{t)hUt) + gfrnhlit)? 

+gf{t)hl{t)hlit) + gUt)[hl{t)]'', (3.45) 

5 

hUt) = Y^airn'r 

r=l 

= gl{t)h\{t)+gl(t)hl{t)+gl(t)[hl{t)? 

+gi{t)hl{t)hl{t) + gl{t)[hl{t)]\ etc. (3.46) 
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And, from (3.7), we have the initial conditions 

K(f) = 51 (3.47) 

We see that if we desire only the degree one terms in the expansion (2.8), then it is 
only necessary to integrate the equations (3.37) through (3.40) with the initial conditions 
(3.47). A moment's reflection shows that so doing amounts to integrating the first-degree 
variational equations. We also observe that if we desire only the degree one and degree 
two terms in the expansion (2.8), then it is only necessary to integrate the equations (3.37) 
through (3.46) with the initial conditions (3.47), etc. 



4 Solution of Complete Variational Equations Using Back- 
ward Integration 

There is another method of determining the /i^ that is surprising, ingenious, and in some 
ways superior to that just described. It involves integrating backward in time [5]. 

4.1 Method of Backward Integration 

Let us rewrite (2.9) in the slightly more explicit form 

Ci = Y.^'a{t\tf)Gr{C) (4.1) 
r 

to indicate that there are two times involved, and . From this perspective, (3.6) is a 
set of differential equations for the quantities (d/dt)h'''^(t'\t) that is to be integrated and 
evaluated at t = t-^. An alternate procedure is to seek a set of differential equations for the 
quantities {d/dt)h^{t,tf) that is to be integrated and evaluated at t = f. 
As a first step in considering this alternative, rewrite (4.1) in the form 

r 

Now reason as follows: If t is varied and at the same time the quantities C(t) are varied 
(evolve) so as to remain on the solution to (2.7) having final conditions then the 
quantities i^^ must remain unchanged. Consequently, there is the differential equation 
result 

= dCVdt = i;[(a/5tl/la(i,tO]G.(C(t)) +E^a(*'*^)(^M")G'r(C(i))- (4.3) 
r r 

Let US introduce the notation }f^{t,t^) for {d/dt)h^{t,tf) so that the first term on the 
right side of (4.3) can be rewritten in the form 

J2[{d/di)hl{t,tf)]Gr{0 = Y,KGr{0. (4.4) 

r r 
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Next, begin working on the second term on the right side of (4.3) by replacing the sum- 
mation index r by the dummy index r', 

j2K{t,tf){d/dt}Grm) = ^h:;it,tf){d/dt)GAat))- (4.5) 

r r' 

Now carry out the indicated differentiation using the chain rule and the relation (2.7) which 
describes how the quantities ( vary along a solution, 

id/dF}Gr'{m) = J2(^Gr'/dCb){dCb/dt) = J2i9Gr'/dCb)gl"it)Gr"{at))- (4-6) 

b br" 

Watch closely: Since the Gr are simply standard monomials in the (, there must be 
relations of the form 

[{d/da)Gr'{0]Gr"{0 = Yl Cl'r"Gr{0 (4-7) 

r 

where the CJ^,^,, are universal constant coefficients that describe certain combinatorial 
properties of monomials. As a result, the second term on the right side of (4.3) can be 
written in the form 

Y,h':{t,tf){d/dt)GAat)) = J2Gr{0 E Cl,,"9t'm:{t,tf). (4.8) 

r' r br'r" 

Since the monomials Gr are linearly independent, the relations (4.3) through (4.8) imply 
the result 

Kit,t^) = - E Cl,r"9b"ma{i,tf)- (4.9) 

br'r" 

This result is a set of differential equations for the that are to be integrated from 
t = t^ hack ioi= f. Also, evaluating (4.2) for i = t^ gives the results 

r 

from which it follows that (with the usual polynomial labeling) the satisfy the final 
conditions 

h:{tf,tf) = 6:. (4.11) 

Therefore the solution to (4.9) is uniquely defined. Finally, it is evident from the definition 
(4.7) that the coefficients Cj^^,^„ satisfy the relation 

CfeV,." = unless [D{r') - 1] + D{r") = D{r). (4.12) 

Therefore, since D{r") > 1 in (4.9), it follows from (4.12) that the only K^^ that occur on 
the right side of (4.9) are those that satisfy 

D{r') < D{r). (4.13) 
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Similarly, the only that occur are those that satisfy 

D{r") < D{r). (4.14) 

Therefore, as before, to determine the coefficients /ij^ of the Taylor expansion (2.9) through 
terms of some degree D, it is only necessary to integrate a finite number of equations, and 
the right sides of these equations involve only the coefficients for this degree and lower. 

Comparison of the differential equation sets (3.6) and (4.9) shows that the latter has 
the remarkable property of being linear in the unknown quantities h'^. This feature means 
that the evaluation of the right side of (4.9) involves only the retrieval of certain universal 
constants Cl^,^,, and straight-forward multiplication and addition. By contrast, the use of 
(3.6) requires evaluation of the fairly complicated nonZmear functions {h^)- Finally, it 
is easier to insure that a numerical integration procedure is working properly for a set of 
linear differential equations than it is for a nonlinear set. 

The only complication in the use of (4.9) is that the equations must be integrated 
backwards in t. Correspondingly the equations (2.6) for the design solution must also be 
integrated backwards since they supply the required quantities (7^ through use of (2.4) and 
(2.5). This is no problem if the final quantities z'^(t^°) are known. However if only the 
initial quantities z'^{t^^) are known, then the equations (2.6) for z*^ must first be integrated 
forward in time to find the final quantities z'^(i^°). 

4.2 The Two-variable Case Revisited 

For clarity, let us also apply this second method to the two- variable case. Table 2 shows 
the nonzero values of C^^,^,, for r G [1,5] obtained using (3.10) through (3.14) and (4.7). 
Note that the rules (4.12) hold. Use of this Table shows that in the two- variable case the 



equations (4.9) take the form 

h\{i,tf) = - g\{t)h\{i,tf)-g\{^hl{i,tf), (4.15) 

hl{i, tf) = - glimW, tf) - gl{F)hlit, tf), (4.16) 

hlit,tf) = - gl{t)h\it,tf)-glit)hl{t,tf), (4.17) 

hUt, tf) = - gfmlit, tf) - gl{t)hl{i, t^), (4.18) 

hUt,tf) = - 2g\{t)h\{ltf) - gl{t)h\{ltf) - g\{f)h\{l^^^ (4.19) 

hl{i, t^) = - 2g\{t)hl{t, tf) - gUt)hl{t, tf) - gl{i)hi{i, t^) - gl{i)hl{i, t^), (4.20) 

-2gl{t)hUt,tf) - gl{t)h\{i,tf) - gl{t)hl{t,tf), (4.21) 
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hiit,tf) = - glit)hl{t,tf)-2gl{t}hUt,tf)-gt{t}hlit,tf) 

-2gUt)hUt,tf) - gl{t)hi{t,tf) - gl{t)hl{t,tf), (4.22) 

hl{t,tf) = - glit)hi{i,tf) - g!it)h\it,tf) - 2gl{t)hl{i,tf) - gl{t)hlit,tf), (4.23) 

hUlt^) = - gfit)hiit,t^) - glmlilff) 2gl{t)hl{i,tf) gUt)hlit,f/), etc. (4.24) 

As advertised, the right sides of (4.15) through (4.24) are indeed simpler than those of 
(3.37) through (3.46). 

Table 2: Nonzero values of Cl^,^„ for r G [1,5] in the two- variable case. 

r b r' r" C 

~\ i i i \~ 

12 2 11 
2 112 1 

2 2 2 2 1 

3 113 1 
3 13 12 
3 2 2 3 1 

3 2 4 1 1 

4 114 1 
4 13 2 2 
4 14 11 
4 2 2 4 1 
4 2 4 2 1 

4 2 5 1 2 

5 115 1 
5 14 2 1 
5 2 2 5 1 
5 2 5 2 1 



5 DufRng Equation Example 
5.1 Introduction 

As an example application, this section studies some aspects of the Duffing equation. 
The behavior of the driven DufHng oscillator, like that of generic nonlinear systems, is 
enormously complicated. Consequently, we will be able to only touch on some of the 
highlights of this fascinating problem. 
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Buffing's equation describes tlie beliavior of a periodically driven damped nonlinear 
oscillator governed by the equation of motion 

X + ax + bx + cx^ = dcos{Clt + ip). (5-1) 

Here ip is an arbitrary phase factor that is often set to zero. For our purposes it is more 
convenient to set 

ip = Tr/2. (5.2) 

Evidently any particular choice of ip simply results in a shift of the origin in time, and this 
shift has no physical consequence since the left side of (5.1) is independent of time. 

We assume b,c > 0, which is the case of a positive hard spring restoring force. ^ We 
make these assumptions because we want the Duffing oscillator to behave like an ordinary 
harmonic oscillator when the amplitude is small, and we want the motion to be bounded 
away from infinity when the amplitude is large. Then, by a suitable choice of time and 
length scales that introduces new variables q and r, the equation of motion can be brought 
to the form 

q + 2/3q + q + q^ = -esiuLOT, (5.3) 

where now a dot denotes d/dr and we have made use of (5.2). In this form it is evident 
that there are 3 free parameters: (3, e, and u. 



5.2 Stroboscopic Map 

While the Duffing equation is nonlinear, it does have the simplifying feature that the 
driving force is periodic with period 

T = 2tt/u}. (5.4) 
Let us convert (5.3) into a pair of first-order equations by making the definition 

p = q, (5.5) 

with the result 

Q = P, 

p = —2f3p — q — q^ — esinwr. (5-6) 

Let q^ ,p^ denote initial conditions at r = 0, and let q'^^p^ be the final conditions resulting 
from integrating the pair (5.6) one full period to the time r = T. Let M denote the transfer 
map that relates q'^TP^ to q^,p^. Then, using the notation z = {q,p), we may write 

z^=Mz^. (5.7) 
^Other authors consider other cases, particularly the 'double well' case & < and c > 0. 
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Suppose we now integrate for a second full period to find q^,p'^- Since the right side of 
(5.6) is periodic, the rules for integrating from r = T to r = 2T are the same as the rules 
for integrating from r = to r = T. Therefore we may write 

z"^ = Mz^ = M^z^, (5.8) 

and in general 

z'^+i =Xz" = X"+iz°. (5.9) 

We may regard the quantities z^ as the result of viewing the motion in the light provided 
by a stroboscope that flashes at the times'^ 

= nT. (5.10) 

Because of the periodicity of the right side of the equations of motion, the rule for sending 
to z""*"^ over the intervals between successive flashes is always the same, namely Ad. For 
these reasons Ai is called a stroboscopic map. Despite the explicit time dependence in the 
equations of motion, because of periodicity we have been able to describe the long-term 
motion by the repeated application of a single fixed map. 

5.3 Feigenbaum Diagram Overview 

One way to study a map and analyze its properties, in this case the Duffing stroboscopic 
map, is to find its fixed points. When these fixed points are found, one can then display how 
they appear, move, and vanish as various parameters are varied. Such a display is often 
called a Feigenbaum diagram. This subsection will present selected Feigenbaum diagrams, 
including the infinite period doubling cascade and its associated strange attractor, for 
the stroboscopic map obtained by high-accuracy numerical integration of the equations of 
motion (5.6). They will be made by observing the behavior of fixed points as the driving 
frequency lo is varied. For simplicity, the damping parameter will be held constant at the 
value P = 0.10. Various sample values will be used for the driving strength e.^ 

Let us begin with the case of small driving strength. When the driving strength is 
small, we know from energy considerations that the steady-state response will be small, 
and therefore the behavior of the steady-state solution will be much like that of the driven 
damped /inear harmonic oscillator. That is, for small amplitude motion, the term in (5.3) 
will be negligible compared to the other terms. We also know that, because of damping, 
there will be only one steady-state solution, and therefore A4 has only one fixed point z-^ 
such that 

Mz^ = zL (5.11) 

^Note that, with the choice (5.2) for ip, the driving term described by the right side of (5.3) vanishes at 
the stroboscopic times r". 

^Of course, one can also make Feigenbaum diagrams in which some other parameter, say e, is varied 
while the others, including cj, are held fixed. 
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Finally, again because of damping, we know that this fixed point is stable. That is, if Ai 
is applied repeatedly to a point near , the result is a sequence of points that approach 
ever closer to . For this reason a stable fixed point is also called an attractor. 

5.3.1 A Simple Feigenbaum Diagram 

Figure 2 shows the values of qf{io) and for the case e = 0.150, and Figure 3 shows pf{uj). In 
the figures the phase-space axes are labeled as goo and poo to indicate that what are being 
displayed are steady-state values reached after a large number of applications of A^. As 
anticipated, we observe from Figures 2 and 3 that the response is much like the resonance 
response of the driven damped linear harmonic oscillator.^ Note that the coefficient of q 
in (5.12) is 1, and therefore at small amplitudes, where can be neglected, the Duffing 
oscillator described by (5.3) has a natural frequency near 1. Correspondingly, Figure 2 
displays a large response when the driving frequency has the value a; ~ 1. Observe, 
however, that the response, while similar, is not exactly like that of the driven damped 
linear harmonic oscillator. For example, the resonance peak at a; ~ 1 is slightly tipped to 
the right, and there is also a small peak for uj ~ 1/3. 

Our strategy for further exploration will be to increase the value of the driving strength 
e, all the while observing the stroboscopic Duffing map Feigenbaum diagram as a function 
of ui. We hasten to add the disclaimer that the driven Duffing oscillator displays an 
enormously rich behavior that varies widely with the parameter values /3, e, uj, and we 
shall be able to give a brief summary of some of it. Also, for brevity, we shall generally 
only display qf{u}). 

5.3.2 Saddle-Node (Blue-Sky) Bifurcations 

Figure 4 shows the q Feigenbaum diagram for the case of somewhat larger driving strength, 
e = 1.50. For this driving strength the resonance peak, which previously occurred at w ~ 1, 
has shifted to a higher frequency and taken on a more complicated structure. There are 
now also noticeable resonances at lower frequencies, with the most prominent one being at 
UJ ~ 1/2. 

Examination of Figure 4 shows that for a; < 1.5 there is a single stable fixed point 
whose trail is shown in black. Then, as uj is increased, a pair of fixed points is born at 
UJ ~ 1.8.^ One of them is stable. The other, whose trail as uj is varied is shown in red, is 
unstable. That is, if A4 is applied repeatedly to a point near this fixed point, the result is 
a sequence of points that move ever farther away from the fixed point. For this reason an 
unstable fixed point is also called a repellor. 

^It was the desire for goo to exhibit a resonance-like peak as a function of uj that dictated the choice 
(5.2) for ijj. 

^Actually, in the analytic spirit of Poincare, these fixed points also exist for smaller values of u, but are 
then complex. They first become purely real, and therefore physically apparent, when a; ~ 1.8. 
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Figure 3: Feigenbaum diagram showing limiting values 

Poo ^ 

function of uj (when /3 = 0.1 

and e = .15) for the stroboscopic Duffing map. 

This appearance of two fixed points out of nowhere is called a saddle-node bifurcation 
or a blue-sky bifurcation, and the associated Fiegenbaum diagram is then sometimes called 
a bifurcation diagram. The original stable fixed point persists as ui is further increased so 
that over some uj range there are 3 fixed points. Then, as u is further increased, the original 
fixed point and the unstable fixed point move until they meet and annihilate when w ~ 2.6.*^ 
This disappearance is called an inverse saddle-node or inverse blue-sky bifurcation. Finally, 
for still larger uj values there is again only one fixed point, and it is stable. 

The appearance and disappearance of stable-unstable fixed-point pairs, as uj is varied, 
has a striking dynamical consequence. Suppose, for example in the case of Figure 4, that 
the driving frequency uj is below the value u; ~ 1.8 where the saddle- node bifurcation 
occurs. Then there is only one fixed point, and it is attracting. Now suppose uj is slowly 
increased. Then, since the fixed-point solution is attracting, the solution for the slowly 
increasing uj case will remain near this solution. See the upper black trail in Figure 4. This 
"tracking" will continue until uj reaches the value uj ~ 2.6 where the inverse saddle-node 
bifurcation occurs. At this value the fixed point being followed disappears. Consequently, 
since the one remaining fixed point is also an attractor, the solution evolves very quickly 
to that of the remaining fixed point. It happens that the oscillation amplitude associated 

^Strictly speaking, a Feigenbaum diagram displays only the trails of stable fixed points while a bifurcation 
diagram displays the trails of all fixed points. 

^Actually, they are not destroyed, but instead become complex and therefore disappear from view. 
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Figure 4: Feigenbaum/bifurcation diagram showing limiting values Qoo as a function of co 
(when /3 = 0.1 and e = 1.5) for the stroboscopic Duffing map. Trails of the stable fixed 
points are shown in black. Also shown, in red, is the trail of the unstable fixed point. 
Finally, jumps in the steady-state amplitude are illustrated by vertical dashed lines at 
ui ~ 1.8 and ui ~ 2.6. 
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with this fixed point is much smaller, and therefore there appears to be a sudden jump 
in oscillation amplitude to a smaller value. Now suppose uj is slowly decreased from a 
value above the value bj ~ 2.6 where the inverse saddle-node bifurcation occurs. Then 
the solution will remain near that of the fixed point lying on the bottom black trail in 
Figure 4. This tracking will continue until a; reaches the value a; ~ 1.8 where the fixed 
point being followed disappears. Again, since the remaining fixed point is attracting, the 
solution will now evolve to that of the remaining fixed point. The result is a jump to 
a larger oscillation amplitude. Evidently the steady-state oscillation amplitude exhibits 
hysteresis as lj is slowly varied back and forth over an interval that begins below the value 
where the first saddle-node bifurcation occurs and ends at a value above that where the 
inverse saddle-node bifurcation occurs. 

5.3.3 Pitchfork Bifurcations 

Let us continue to increase e. Figure 5 shows that a qualitatively new feature appears 
when e is near 2.2: a bubble is formed between the major resonant peak (the one that 
has saddle-node bifurcated) and the subresonant peak immediately to its left. To explore 
the nature of this bubble, let us make e still larger, which, we anticipate, will result in 
the bubble becoming larger. Figure 6 shows the Feigenbaum diagram in the case e = 5.5. 
Now the major resonant peak and the subresonant peak have moved to larger u values. 
Correspondingly, the bubble between them has also moved to larger oj values. Moreover, 
it is larger, yet another smaller bubble has formed, and the subresonant peak between 
them has also undergone a saddle-node bifurcation. For future use, we will call the major 
resonant peak the first or leading saddle-node bifurcation, and we will call the subresonant 
peak between the two bubbles the second saddle-node bifurcation, etc. Also, we will call 
the bubble just to the left of the first saddle- node bifurcation the first or leading bubble, 
and the next bubble will be called the second bubble, etc. 

We also note that three short trails have appeared in Figure 6 just to the right of 
OJ = A. They correspond to a period-three bifurcation followed shortly thereafter by an 
inverse bifurcation. Actually, much closer examination shows that there are six trails 
consisting of three closely-spaced pairs. Each pair comprises a stable and an unstable fixed 
point of the map M^. They are not fixed points of M. itself, but rather are sent into each 
other in cyclic fashion under the action of M. 

Figure 7 shows the larger (leading) bubble in Figure 6 in more detail and with the 
addition of red lines indicating the trails of unstable fixed points. It reveals that the bubble 
describes the simultaneous bifurcation of a single fixed point into three fixed points. Two 
of these fixed points are stable and the third, whose q coordinate as a function of cu are 
shown as a red line, is unstable. What happens is that, as uj is increased, a single stable 
fixed point becomes a triplet of fixed points, two of which are stable and one of which is 
unstable. This is called a pitchfork bifurcation. Then, as oj is further increased, these three 
fixed points again merge, in an inverse pitchfork bifurcation, to form what is again a single 
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Figure 6: Feigenbaum diagram showing limiting values goo as a function of oo (when /3 = 0.1 
and e = 5.5) for the stroboscopic Duffing map. The first bubble has grown, a second smaller 
bubble has formed to its left, and the sub-resonant peak between them has saddle- node 
bifurcated to become the second saddle-node bifurcation. 
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5.3.4 A Plethora of Bifurcations and Period Doubling Cascades 

We end our study of the Duffing equation by increasing e from its earlier value e = 5.5 to 
much larger values. First we will set e = 22.125. Based on our experience so far, we might 
anticipate that the Feigcnbaum diagram would become much more complicated. That is 
indeed the case. Figure 8 displays q^o when /? = 0.1 and e = 22.125, as a function of oj, for 
the range w G (0, 12). Evidently the behavior of the attractors for the stroboscopic Duffing 
map, which is what is shown in Figure 8, is extremely complicated. There are now a great 
many fixed points both of Ai itself and various powers of M. 

Of particular interest to us are the two areas around w = .8 and oj = 1.25. They 
contain what has become of the first two bubbles in Figure 7, and are shown in greater 
magnification in Figure 9. What has happened is that bubbles have formed within bubbles, 
and bubbles have formed within these bubbles, etc. to form a cascade. However, these 
interior bubbles are not the result of pitchfork bifurcations, but rather the result of period- 
doubling bifurcations. For example, the bifurcation that creates the first bubble at a; ~ 1.2 
is a pitchfork bifurcation. But the successive bifurcations within the bubble are period- 
doubling bifurcations. In a period-doubling bifurcation a fixed point that is initially stable 
becomes unstable as lj is increased. When this happens, simultaneously two stable fixed 
points of Ai^ are born. They are not fixed points of Ai itself, but rather are sent into 
each other under the action of M. Hence the name "period doubling". The map M. must 
be applied twice to send such a fixed point back into itself. In the next period doubling, 
fixed points of A4'^ are born, etc. However we note that, as to increases, the sequence of 
period-doubling bifurcations only occurs a finite number of times and then undoes itself. 

Remarkably, when e is just somewhat larger, infinite sequences of period doubling 
cascades can occur. Figure 10 shows what happens when e = 25. 

5.3.5 More Detailed Vievi^ of Infinite Period Doubling Cascades 

To display the infinite period doubling cascades in more detail. Figure 11 shows an en- 
largement of part of Figure 10. And Figures 12 and 13 show successive enlargements of 
parts of the first bubble in Figure 10. From Figure 11 we see that the first bubble forms as 
a result of a pitchfork bifurcation just to the right of cj = 1.2, and from Figures 12 and 13 
we see that the first period doubling bifurcation occurs in the vicinity oi u = 1.268. From 
Figure 13 it is evident that successive period-doublings occur an infinite number of times 
to ultimately produce a chaotic region when ui exceeds uj 2± 1.29. 

5.3.6 Strange Attractor 

As evidence that the behavior in this region is chaotic. Figures 14 and 15 show portions 
of the full phase space, the q,p plane, when u = 1.2902. Note the evidence for fractal 
structure. The points appear to lie on a strange attractor. 
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Figure 8: Feigenbaum diagram showing limiting values Qoo as a function of a; (when /3 = 0.1 
and e = 22.125) for the stroboscopic Duffing map. 
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Figure 9: Enlarged portion of the Feigenbaum diagram of Figure 8 displaying limiting 
values goo as a function of u (when f3 = 0.1 and e = 22.125) for the stroboscopic Duffing 
map. It shows part of the first bubble at the far right, the second bubble, and part of a 
third bubble at the far left. Examine the first and second bubbles. Each initially consists 
of two stable period-one fixed points. Each also contains the beginnings of period-doubling 
cascades. These cascades do not complete, but rather remerge to again result in a pair of 
stable period-one fixed points. There are also many higher-period fixed points and their 
associated cascades. 
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Figure 11: Enlargement of a portion Figure 10 showing the first, second, and third bubbles. 
The period-doubling cascades in each of the first and second bubbles complete. Then they 
undo themselves as uj is further increased. There is no period doubling in the third bubble 
when e = 25. 
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Figure 12: Detail of part of the first bubble in Figure 11 showing upper and lower infinite 
period-doubling cascades. Part of the trail of the stable fixed point associated with the 
second saddle-node bifurcation accidentally appears to overlay the upper period doubling 
bifurcation. Finally, there are numerous cascades and remergings associated with higher- 
period fixed points. 
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Figure 13: Detail of part of the upper cascade in Figure 12 showing an infinite period- 
doubling cascade, followed by chaos, for what was initially a stable period-one fixed point. 
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Figure 14: Limiting values of goo, Poo for the stroboscopic Duffing map when 00 = 1.2902 
(and P = .1 and e = 25). They appear to lie on a strange attractor. 
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Figure 15: Enlargement of boxed portion of Figure 14 illustrating the beginning of self- 
similar fractal structure. 
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6 Polynomial Approximation to Duffing Stroboscopic Map 



In this section we will find the complete variational equations for the Duffing equation 
including dependence on the parameter lo. The two remaining parameters, P and e, will 
remain fixed. We will then study the properties of the resulting polynomial approximation 
to the Duffing stroboscopic map obtained by truncating its Taylor expansion. 

6.1 Complete Variational Equations with Selected Parameter Depen- 
dence 

To formulate the complete variational equations we could proceed as in Section 2.2 by 
setting uj = and zs = w'^ + ^3- So doing would require expanding the function sin(a;T) = 
sin[(a;'^ + C3)t] as a power series in Cs- Such an expansion is, of course, possible, but it 
leads to variational equations with a large number of forcing terms g^^ since the expansion 
of sin[(w'^ + Cs)"?"] contains an infinite number of terms. 

In the case of Duffing's equation this complication can be avoided by a further change 
of variables. Recall that the first change of variables brought the Duffing equation to the 
form (5.3), which we now slightly rewrite as 

q" + 2Pq' + q + q^ = -esinuT (6.1) 

where d/dr is now denoted by a prime. Next make the further change of variables 

q = uQ, (6.2) 

oj = 1/(7, (6.3) 
OJT = t. (6.4) 
When this is done, there are the relations 

q' = LJ^Q (6.5) 

and 

q" = uj^Q (6.6) 

where now a dot denotes d/dt. [Note that the variable t here is different from that in (5.1).] 
Correspondingly, Duffing's equation takes the form 

Q + 2l5aQ + a'^Q + Q^ = -ea^ sin t. (6.7) 

This equation can be converted to a first-order set of the form (2.1) by writing 

Q = zx (6.8) 
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and 

Q = Z2. (6.9) 

Following the method of Section 2.2, we augment the first-order equation set associated 
with (6.7) by adding the equation 

cr = 0. (6.10) 

Then we may view a as a variable, and (6.10) guarantees that this variable remains a 
constant. Taken together, (6.7) and (6.10) may be converted to a first-order triplet of the 
form (2.12) by writing (6.8), (6.9), and 

(7 = ^3- (6.11) 

Doing so gives the system 

zi = Z2, (6.12) 
Z2 = —2/3^322 — z'^zi — zf — ezf sint, (6.13) 
is = 0, (6.14) 

and we see that there are the relations 

fiiz,t) = Z2, (6.15) 

f2{z,t) = - 2PZ3Z2 - zjzi -zf- ezlsint, (6.16) 

f3{z,t) = 0. (6.17) 

Note that, with the change of variables just made, the stroboscopic map is obtained by 
integrating the equations of motion from t = to t = 27r. 

As before, we introduce deviation variables using (2.2) and carry out the steps (2.3) 
through (2.9). In particular, we write 

Z3 = zl + Cs = CJ'^ + Cs. (6.18) 

This time we are working with monomials in the three variables <^i, (^2, and (^3. [That is, 
a ranges from 1 to 3 in (2.12).] They are conveniently labeled using the indices r given in 
Table 3 below. With regard to the expansions (2.4) and (2.5), we find the results 

/i(z'^ + C,0 = 4 + C2, (6.19) 



f2{z'' + C,t) = -2/3(zf + C3)(4 + C2)-(4 + C3)'(^f + Cl) 

-{zf + Cif- eizi + (3? sint 
= [-2^zlzi - zfizif - {zff - e{zif sint] 

-[3(zf )2 + (zf )2]Ci - 2(3ziQ2 - Wzi + 2zizi + ?>e{zif sintjCs 
-2/3C2C3 - {zf + 3ezf )C| - C1C2 - 3zfCi 

-C?-ClC|-e(sin^)C3^ (6.20) 
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/3(z'^ + C,t) = 0. (6.21) 

Note the right sides of (6.19) through (6.21) are at most cubic in the deviation variables 
(a- Therefore, from Table 3, we see that the index r for the g'^ should range from 1 through 
19. It follows that for Buffing's equation (with a parameter expansion) the only nonzero 
forcing terms are given by the relations 



5? = 1, 


(6.22) 


gl = -Sizff - (zif, 


(6.23) 


gl = -Wzi 


(6.24) 


-2^4 - '^44 - 3e(zf)2sini, 


(6.25) 


gi = '^zf, 


(6.26) 


gl = -'^zi, 


(6.27) 


gl = -2/3, 


(6.28) 


f2 = -zf-Sezismt, 


(6.29) 


gl' = -1, 


(6.30) 


g¥ = -1, 


(6.31) 


gl^ = — esini. 


(6.32) 
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Table 3: A labeling scheme for monomials in three variables. 
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6.2 Performance of Polynomial Approximation 

Let Als denote the 8*'^-order polynomial map (with parameter dependence) approximation 

to the stroboscopic Duffing map M. Provided the relevant phase-space region is not too 
large, we have found that A^g reproduces all the features, described in Section 5.3, of 
the exact map[6]. (The phase-space region must lie within the convergence domain of the 
Taylor expansion.) This reproduction might not be too surprising in the cases of elementary 
bifurcations such as saddle-node and pitchfork bifurcations. What is more fascinating, as 
we will see, is that M.^ also reproduces the infinite period doubling cascade and associated 
strange attractor. 

6.2.1 Infinite Period Doubling Cascade 

Figure 16 shows the partial Feigenbaum diagram for the map in the case that ^9 = 0.1 
and e = 25. The Taylor expansion is made about the point indicated by the black dot. 
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This point has the coordinates 

gbd = 1.26082, pbd = 2.05452, Wbd = 1-285. (6.33) 

It was selected to be an unstable fixed point of Ai, but that is not essential. Any nearby 
expansion point would have served as well. Note the remarkable resemblance of Figures 13 
and 16. 

We have referred to Figure 16 as a partial Feigenbaum diagram because it shows only 
Qoo and not poo- In order to give a complete picture. Figure 17 displays them both. 




Figure 16: Partial Feigenbaum diagram for the map A4s- The black dot marks the point 
about which jH is expanded to yield Aig 
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Figure 17: Full Feigenbaum diagram for the map TWs- The black dot again marks the 
expansion point. 
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6.2.2 Strange Attractor 



As displayed in Figures 18 through 21 Ais, hke Ai, appears to have a strange attractor. 
Note the remarkable agreement between Figures 14 and 15 for M and their Ms coun- 
terparts, Figures 18 and 19. In the case of A^g we have been able to obtain additional 
enlargements, Figures 20 and 21, further illustrating a self-similar fractal structure. Analo- 
gous figures are more difficult to obtain for the exact map Ai due to the excessive numerical 
integration time required. By contrast the map Ms, because it is a simple polynomial, is 
easy to evaluate repeatedly. 
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Figure 18: Limiting values of QocPoo for the map Ms when cj = 1.2902. They appear to 
lie on a strange attractor. 
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Figure 20: Enlargement of boxed portion of Figure 19 illustrating the continuation of 
self-similar fractal structure. 
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Figure 21: Enlargement of boxed portion of Figure 20 illustrating the further continuation 
of self-similar fractal structure. 
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7 Numerical Implementation 



The forward integration method (Section 3.1) can be implemented by a code employing 
the tools of automatic differentiation (AD) described by Neidinger [7].*^ In this approach 
arrays of Taylor coefficients of various functions are referred to as AD variables or pyramids 
since, as will be seen, they have a hyper- pyramidal structure. Generally the first entry in 
the array will be the value of the function about some expansion point, and the remaining 
entries will be the higher-order Taylor coefficients about the expansion point and truncated 
beyond some specified order. Such truncated Taylor expansions are also commonly called 
jets. 

In our application elements in these arrays will be addressed and manipulated with the 
aid of scalar indices associated with look-up tables generated at run time. We have also 
replaced the original APL implementation of Neidinger with a code written in the language 
of Mathematica (Version 6, or 7) [8]. Where necessary, for those unfamiliar with the details 
of Mathematica, we will explain the consequences of various Mathematica commands. The 
inputs to the code are the right sides (RS) of (1.1), or (2.1). Other input parameters are 
the number of variables m, the desired order of the Taylor map p, and the initial conditions 
{z^y for the design-solution equation (2.3). 

Various AD tools for describing and manipulating pyramids are outlined in Section 
7.1. There we show how pyramid operations are encoded in the case of polynomial RS, 
as needed for the Duffing equation. For brevity, we omit the cases of rational, fractional 
power, and transcendental RS. These cases can also be handled using various methods 
based on functional identities and known Taylor coefficients, or the differential equations 
that such functions obey along with certain recursive relations [7]. In Section 7.2, based on 
the work of Section (7.1), we in effect obtain and integrate numerically the set of differential 
equations (3.6) in pyramid form, i.e. valid for any map order and any number of variables. 
Section 7.3 treats the specific case of the Duffing equation. A final Section 7.4 describes 
in more detail the relation between integrating equations for pyramids and the complete 
variational equations. 

7.1 AD tools 

This section describes how arithmetic expressions representing fa{z,t), the right sides of 
(1.1) where z denotes the dependent variables, are replaced with expressions for arrays 
(pyramids) of Taylor coefficients. These pyramids in turn constitute the input to our code. 
Such an ad- hoc replacement, according to the problem at hand, as opposed to operator 
overloading where the kind of operation depends on the type of its argument, is also the 
approach taken in [7,9,10]. 

^Some authors refer to AD as truncated power series algebra (TPSA) since AD algorithms arise from 
manipulating multivariable truncated power series. Other authors refer to AD as Differential Algebra (DA). 
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Let u,v,w be general arithmetic expressions, i.e. scalar-valued functions of z. They 
contain various arithmetic operations such as addition, multiplication (*), and raising to a 
power (a). (They may also entail the computation of various transcendental functions such 
as the sine function, etc. However, as stated earlier, for simplicity we will omit these cases.) 
The arguments of these operations may be a constant, a single variable or multiple variables 
Za, or even some other expression. The idea of AD is to redefine the arithmetic operations 
in such a way (sec Definition 1), that all functions u, v, w can be consistently replaced with 
the arrays of coefficients of their Taylor expansions. For example, by redefining the usual 
product of numbers (*) and introducing the pyramid operation PROD, u*v is replaced with 
PROD [U,V]. 

We use upper typewriter font for pyramids (U,V, . . .) and for operations on pyramids 
(PROD, POW, ...). Everywhere, equalities written in typewriter fonts have equivalent 
Mathematica expressions. That is, they have associated realizations in Mathematica and 
directly correspond to various operations and commands in Mathematica. In effect, our 
code operates entirely on pyramids. However, as we will see, any pyramid expression 
contains, as its first entry, its usual arithmetic counterpart. 

We begin with a description of our method of monomial labeling. In brief, we list 
all monomials in a polynomial in some sequence, and label them by where they occur in 
the list. Next follow Definition 1 and the recipes for encoding operations on pyramids. 
Subsequently, by using Definition 2, which simply states the rule by which an arithmetic 
expressions is replaced with its pyramid counterpart, we show how a general expression 
can be encoded by using only the pyramid of a constant and of a single variable. 

7.1.1 Labeling Scheme 

A monomial Gj{z) in m variables is of the form 

Gj{z) = {z,y^{z2y^---{zmy-. (7.1) 

Here we have introduced an exponent vector j by the rule 

3 = {ji,j2,---jm)- (7.2) 

Evidently j is an m-tuple of non-negative integers. The degree of Gj(z), denoted by 
is given by the sum of exponents, 

|j|=Jl+j2H him- (7.3) 

The set of all exponents for monomials in m variables with degree less than or equal to p 
will be denoted by r^, 

r?n = {j\\j\<p}. (7.4) 
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It can be shown that this set has L(m,p) entries with L(m,p) given given by a binomial 
coefficient, 

L(™,p) =(";"'). (7.5) 

In [6] this quantity is called 5o(m,p). Assuming that m and p are fixed input variables, we 
will often write F and L. With this notation, a Taylor series expansion (about the origin) 
of a scalar- valued function u of m variables z = {zi,Z2-i ■ ■ ■ Zm), truncated beyond terms of 
degree p, can be written in the form 

u{z) = J2 Gj{z). (7.6) 

Here, for now, U simply denotes an array of numerical coefficients. When employed in code 
that has symbolic manipulation capabilities, each U(jf) may also be a symbolic quantity. 

To proceed, what is needed is some way of listing monomials systematically. With 
such a list, as already mentioned, we may assign a label r to each monomial based on 
where it appears in the list. A summary of labeling methods, and an analysis of storage 
requirements, may be found in [6,9]. Here we describe one of them that is particularly 
useful for our purposes. 

The first step is to order the monomials or, equivalently, the exponent vectors. One 
possibility is lexicographic {lex) order. Consider two exponent vectors j and k. Let j — k 
be the vector whose entries are obtained by component-wise subtraction, 

j -k = {ji - ki,j2 -k2,--- jm - km)- (7.7) 

We say that the exponent vector j is lexicographically greater than the exponent vector 
k, and write j >iex k, if the left-most nonzero entry in j — fc is positive. Thus, for example 
in the case of monomials in three variables {zi,Z2,Z3) with exponents (^1,^2, is), we have 
the ordering 

(1,0,0) >iex (0,1,0) >iex (0,0,1) (7.8) 

and 

(4,2,l)>iex(4,2,0)>iex(2,5,l). (7.9) 

For our purposes we have found it convenient to label the monomials in such a way that 
monomials of a given degree D = \j\ occur together. One possibility is to employ graded 
lexicographic (glex) ordering. If j and k are two exponent vectors, we say that j >giex k if 
either \j\ > \k\, or \j\ = \k\ and j >iex k. 

Table 4 shows a list of monomials in three variables. As one goes down the list, first 
the monomial of degree D = appears, then the monomials of degree D = 1, etc. Within 
each group of monomials of fixed degree the individual monomials appear in descending 
lex order. Note that Table 4 is similar to Table 3 except that it begins with the monomial 
of degree 0. 
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Table 4: A labeling scheme for monomials in three variables. 
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We give the name modified glex sequencing to a monomial listing of the kind shown in 
Table 4. This is the labeling scheme we will use. Other possible listings include ascending 
true glex order in which monomials appear in ascending lex order within each group of 
degree D, and lex order for the whole monomial list as in [7]. 

With the aid of the scalar index r the relation (7.6) can be rewritten in the form 

L{m,p) 

u{z) = U(r)G,(z), (7.10) 

r=l 
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because (by construction and with fixed m) for each positive integer r there is a unique 
exponent j{r), and for each j there is a unique r. Here U may be viewed as a vector with 
entries U(r), and Gr{z) denotes Gj(^r)i^)- 

Consider, in an m-dimensional space, the points defined by the vectors j £ Tm- See 
(7.4). Figure 22 displays them in the case m = 3 and p = 4. Evidently they form a grid 
that lies on the surface and interior of what can be viewed as an m-dimensional pyramid in 
m-dimensional space. At each grid point there is an associated coefficient U(r). Because of 
its association with this pyramidal structure, we will refer to the entire set of coefficients 
in (7.6) or (7.10) as the pyramid U of u{z). 
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Figure 22: A grid of points representing the set Fg. For future reference a subset of Fg, 
called a hox, is shown in blue. 



7.1.2 Implementation of Labeling Scheme 

We have seen that use of modified glex sequencing, for any specified number of variables 
m, provides a labeling rule such that for each positive integer r there is a unique exponent 
j(r), and for each j there is a unique r. That is, there is a invertible function r{j) that 
provides a 1-to-l correspondence between the positive integers and the exponent vectors j. 
To proceed further, it would be useful to have this function and its inverse in more explicit 
form. 

First, we will learn that there is a formula for r(jf), which we will call the Giorgilli 
formula [6]. For any specified m the exponent vectors j take the form (7.2) where all the 
entries ji are positive integers or zero. Begin by defining the integers 

i-\ 

Jm) = i- 1 + (7.11) 
k=0 
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for i G {l,2,---m}. Then, to the general monomial Gj{z) or exponent vector j, assign 
the label 

m 

r{j) = r{ji, ■ . ■ j^) = 1 + Binomial [n(£;ii, • --jm),^]- (7.12) 

e=i 

Here the quantities 

Binomial [n,e]=( ) = { ' ^ ^ ^ ^ ^ I (7.13) 

^ ' ^ V ^ / \ , otherwise J ^ ^ 

denote the usual binomial coefficients. It can be shown that this formula reproduces the 
results of modified glex sequencing [6] . 

Below is simple Mathematica code that implements the Giorgilli formula in the case of 
three variables, and evaluates it for selected exponents j. Observe that these evaluations 
agree with results in Table 4. 

Gfor[jl_,j2_,j3_] := ( 

sl = j3; s2 = l + j3 + j2; s3 = 2 + jS + j2 + j 1; 

tl = Binoniial[sl, l]; t2 = Binoinial[s2, 2]; tS = Binoinial[s3, 3]; 

r = 1 + tl + t2 + t3; r 

) 

Gf or[0, 0, 0] 

Gfor[l,0,0] 

Gfor[2,0, 1] 

Gfor[l,2, 1] 

1 

2 

13 

28 (7.14) 

Second, for the inverse relation, we have found it convenient to introduce a rectangular 
matrix associated with the set Fm. By abuse of notation, it will also be called F. It has 
L{m,p) rows and m columns with entries 

rr,a=ja{r). (7.15) 

For example, looking a Table 4, we see (when m = 3) that Fi^i = and Fi7_2 = 3. Indeed, 
if the first and last columns of Table 4 are removed, what remains (when m = 3) is the 
matrix F^^a- In the language of [6], F is a look up table that, given r, produces the associated 
j. In our Mathematica implementation F is the matrix GAMMA with elements GAMMA[[r, a]]. 
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The matrix GAMMA is constructed using the Mathematica code illustrated below, 

Needs [" Comb inat or ica'"]; 
m = 3; p = 4; 

GAMMA = Compositions [0,in]; 

Do[GAMMA = Join[GAMMA, Reverse [Compositions [d,m]]], {d, l,p, l}]; 

L = Length[GAMMA] 

r = 17;a = 2; 

GAMMA[[r]] 

GAMMA[[r, a]] 

35 

{0,3,0} 

3 (7.16) 

It employs the Mathematica commands Compositions, Reverse, and Join. 

We will first describe the ingredients of this code and illustrate the function of each: 

• The command Needs["Combinatorica"']; loads a combinatorial package. 

• The command Compositions[i, m] produces, as a list of arrays (a rectangular array), 
all compositions (under addition) of the integer i into m integer parts. Further- 
more, the compositions appear in ascending lex order. For example, the command 
Compositions[0, 3] produces the single row 

(7.17) 

As a second example, the command Compositions[l, 3] produces the rectangular 
array 

1 

1 

1 (7.18) 

As a third example, the command Compositions[2, 3] produces the rectangular array 

2 

1 1 

2 

1 1 

1 1 

2 (7.19) 
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• The command Reverse acts on the hst of arrays, and reverses the order of the Hst 
while leaving the arrays intact. For example, the nested sequence of commands 
Reverse[Compositions[l, 3]] produces the rectangular array 

1 

1 

1 (7.20) 

As a second example, the nested sequence of commands Reverse[Compositions[2, 3]] 
produces the rectangular array 

2 

1 1 
1 1 
2 
1 1 

2 (7.21) 

Now the compositions appear in descending lex order. 

• Look, for example, at Table 4. We sec that the exponents ja for the r = 1 entry 
are those appearing in (7.17). Next, exponents for the r = 2 through r = 4 entries 
are those appearing in (7.20). Following them, the exponents for the r = 5 through 
r = 10 entries, are those appearing in (7.21), etc. Evidently, to produce the exponent 
list of Table 4, what we must do is successively join various lists. That is what the 
Mathematica command Join can accomplish. 

We are now ready to describe how GAMMA is constructed: 

• The second line in (7.16) sets the values of m and p. They are assigned the values 
m = 3 and p = 4 for this example, which will construct GAMMA for the case of Table 4. 
The third line in (7.16) initially sets GAMMA to a row of m zeroes. The fourth line is 
a Do loop that successively redefines GAMMA by generating and joining to it successive 
descending lex order compositions. The net result is the exponent list of Table 4. 

• The quantity L = L(m, p) is obtained by applying the Mathematica command Length 

to the the rectangular array GAMMA. 

• The last 6 lines of (7.16) illustrate that L is computed properly and that the com- 
mand GAMMA[[r, a]] accesses the array GAMMA in the desired fashion. Specifically, in 
this example, we find from (7.5) that L{3,4) = 35 in agreement with the Mathe- 
matica output for L. Moreover, GAMMA[[17]] produces the exponent array {0,3,0}, in 
agreement with the r = 17 entry in Table 4, and GAMMA[[17, 2]] produces Ti7^2 = 3, 
as expected. 
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7.1.3 Pyramid Operations: Addition and Multiplication 

Here we derive the pyramid operations in terms of j- vectors by using the ordering previously 
described, and provide scripts to encode them in the r-representation (7.10). 

Definition 1 Suppose thatw{z) arises from carrying out various arithmetic operations on 
u{z) and v{z), and the pyramids U and V are known. The corresponding pyramid operation 
on U and V is so defined that it yields the pyramid W of w{z). 

Here we assume that u,v,w are polynomials such as (7.6). 

Wc begin with the operations of scalar multiplication and addition, which are easy to 
implement. If 

w{z) = cu{z), (7.22) 

then 

W(r) = cU(r), (7.23) 

and we write 

W = c U. (7.24) 

If 

w{z) = u{z) +v{z), (7.25) 

then 

W(r) =U(r) +V(r), (7.26) 

and we write 

W = U + V. (7.27) 

In both cases all operations are performed coordinate- wise (as for vectors). 

Implementation of scalar multiplication and addition is easy in Mathematica because, as 
the example below illustrates, it has built in vector routines. There we define two vectors, 
multiply them by scalars, and add the resulting vectors. 

Unprotect[V]; 
U = {1,2,3}; 

V = {4,5,6}; 
W= .1U+ .2V 

{.9,1.2,1.5} (7.28) 

Since V is a "protected" symbol in the Mathematica language, and, for purposes of illus- 
tration, we wish to use it as an ordinary vector variable, it must first be unprotected as in 
line 1 above. The last line shows that the Mathematica output is indeed the desired result. 
The operation of polynomial multiplication is more involved. Now we have the relation 

w{z) =u{z)*v{z), (7.29) 
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and we want to encode 

W = PROD[U,V]. (7.30) 

Let us write u{x) in the form (7.6), but with a change of dummy indices, so that it has 
the representation 

u{z) = U(i) Giiz). (7.31) 



Similarly, write v(z) in the form 



v{z) = J2 Gj{z). (7.32) 



Then, according to Leibniz, there is the result 

u{z)*viz)= E U(i)V(j)Gi(z)*G,-(z). (7.33) 

From (7.1) we observe that 

Gi{z) * Gjiz) = {zir {Z2Y' ■ ■ ■ {zmY- * (ziy^ {Z2Y' ■ ■ ■ {z„,y-^ 

= {ziy^+^\z2r^^' ■ ■ ■ izmy--^'- = Gi+j{z). (7.34) 
Therefore, we may also write 

uiz)*viz)= Y E ^iinj)Gi+j{z). (7.35) 

Now we see that there are two complications. First, there may be terms on the right side 
of (7.35) whose degree is higher than p and therefore need not be computed. Second, there 
are generally many terms on the right side of (7.35) that contribute to a given monomial 
term in w{z) = u{z) * v{z). Suppose we write 

w{z) = Y^{k)Gk{z). (7.36) 

k 

Upon comparing (7.35) and (7.36) we conclude that 

W(fc) = Y U(i)V(j) = E U(fe - mJ)- (7.37) 

i+j=k j<k 

Here, by j < k, we mean that the sum ranges over all j such that ja < ka for all a G [1, m]. 
That is, 

j < k <^ ja < ka for all a G [1, m]. (7.38) 
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Evidently, to implement the relation (7.37) in terms of r labels, we need to describe 
the exponent relation j < k in terms of r labels. Suppose k is some exponent vector with 
label r{k) as, for example, in Table 4. Introduce the notation 

k = r{k). (7.39) 

This notation may be somewhat confusing because k is not the norm of the vector k, but 
rather the label associated with k. However, this notation is very convenient. Now, given 
a label k, we can find k. Indeed, from (7.15), we have the result 

ka = rk,a. (7.40) 

Having found k, we define a set of exponents by the rule 

Bk = {j\j < k}. (7.41) 

This set of exponents is called the k^^ box. For example (when m = 3), suppose k = 28. 
Then we see from Table 4 that k = (1, 2, 1). Table 5 lists, in modified glex order, all the 
vectors in B28, i-e. all vectors j such that j < (1,2,1). These are the vectors shown in 
blue in Figure 22. Finally, with this notation, we can rewrite (7.37) in the form 

W(/c)= ^ U(fe-j)V(j). (7.42) 



Table 5: The vectors in S28 = {j |i < (1, 2, 1)}. 
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What can be said about the vectors {k — j) as j ranges over Bg? Table 6 lists, for 
example, the vectors j G B28 and the associated vectors i with i = {k — j). Also listed are 
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the labels r(j) and r(i). Compare columns 2,3,4, which specify the j G ^28, with columns 
5,6,7, which specify the associated i vectors. We see that every vector that appears in the 
j list also occurs somewhere in the i list, and vice versa. This to be expected because the 
operation of multiplication is commutative: we can also write (7.37) in the form 

W(fc)= ^ U(j)V(fc-j). (7.43) 

We also observe the more remarkable feature that the two lists arc reverses of each other: 
running down the j list gives the same vectors as running up the i list, and vice versa. 
This feature is a consequence of our ordering procedure. 

As indicated earlier, what we really want is a version of (7.37) that involves labels 
instead of exponent vectors. Looking at Table 6, we see that this is easily done. We may 
equally well think of Bk as containing a collection of labels r{j), and we may introduce a 
reversed array Brevk of complementary labels r'^{j) where 

r'{j) = r{t). (7.44) 

That is, for example, B28 would consist of the first column of Table 6 and Breves would 
consist of the last column of Table 6. Finally, we have already introduced k as being the 
label associated with k. We these understandings in mind, we may rewrite (7.37) in the 
label form 

W(A;) = ^ U(r'^)V(r) = ^ U(r)V(r^). (7.45) 

This is the rule W = PROD[U, V] for multiplying pyramids. In the language of [6], B^ and 
Brevk are look hack tables that, given a k, look back to find all monomial pairs with labels 
r, t"^ which produce, when multiplied, the monomial with label k. 

7.1.4 Implementation of Multiplication 

The code shown below in (7.46) illustrates how Bk and Brevk are constructed using Math- 
ematica. 

JSK[list_,K_] := 

Position[Apply[And, Thread[#l<=#2& [#,K] ] ] & /@ list, True] //Flatten; 
B = Table[JSK[GAMMA, GAmA[[k]]], {k, 1, L}]; 

Brev = Reverse /@ B; (7.46) 

As before, some explanation is required. The main tasks are to implement the j < k 
operation (7.38) and then to employ this implementation. We will begin by implementing 
the j < k operation. Several steps are required, and each of them is described briefly 
below: 
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Table 6: The vectors j and i = {k — j) for j G B28 and ka = ^28,a- 
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• When Mathematica is presented with a statement of the form j <= k, with j and 
k being integers, it repHes with the answer True or the answer False. (Here j <= k 
denotes j < k.) Two sample Mathematica runs are shown below: 

3 <= 4 

TVue (7.47) 
5 <= 4 

False (7.48) 

• A Mathematica function can be constructed that does the same thing. It takes the 
form 

#l<=#2&[j,k] (7.49) 

Here the symbols #1 and #2 set up two slots and the symbol & means the operation 

to its left is to be regarded as a function and is to be applied to the arguments to 
its right by inserting the arguments into the slots. Below is a short Mathematica run 
illustrating this feature. 

j =3;k = 4; 
#1 <= #2 & [j,k] 

True (7.50) 
Observe that the output of this run agrees with that of (7.47). 
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• The same operation can be performed on pairs of arrays (rather than pairs of num- 
bers) in such a way that corresponding entries from each array arc compared, with 
the output then being an array of True and False answers. This is done using the 
Mathematica command Thread. Below is a short Mathematica run illustrating this 
feature. 

j ={l,2,3};k = {4,5,l}; 
Thread[#l <= #2 & [j,k]] 

{True, True, False} (7-51) 

Note that the first two answers in the output array are True because the statements 
1 < 4 and 2 < 5 are true. The last answer in the output array is False because the 
statement 3 < 1 is false. 

• Suppose, now, that we are given two arrays j and k and we want to determine if 
j < k in the sense of (7.38). This can be done by applying the logical And operation 
(using the Mathematica command Apply) to the True/False output array described 
above. Below is a short Mathematica run illustrating this feature. 

j ={l,2,3};k = {4,5,l}; 
Apply[And, Thread[#l <= #2 & [j,k]]] 

False (7.52) 

Note that the output answer is False because at least one of the entries in the output 
array in (7.51) is False. The output answer would be True if, and only if, all entries 
in the output array in (7.51) were True. 

• Now that the j < k operation has been defined for two exponent arrays, we would 
like to construct a related operator/function, to be called JSK. (Here the letter S 
stands for smaller than or equal to.) It will depend on the exponent array k, and its 
task will be to search a list of exponent arrays to find those j within it that satisfy 
j < k. The first step in this direction is to slightly modify the function appearing in 
(7.52). Below is a short Mathematica run that specifies this modified function and 
illustrates that it has the same effect. 

j ={l,2,3};k = {4,5,l}; 

Apply[And,Thread[#l <= #2 & [#,k]]] & [j] 

False (7.53) 

Comparison of the functions in (7.52) and (7.53) reveals that what has been done is to 

replace the argument j in (7.52) by a slot #, then follow the function by the character 
&, and finally add the symbols [j]. What this modification does is to redefine the 
function in such a way that it acts on what follows the second &. 
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• The next step is to extend the function appearing in (7.53) so that it acts on a hst of 
exponent arrays. To do this, we replace the symbols [j] by the symbols /@ list. The 
symbols /@ indicate that what stands to their left is to act on what stands to their 
right, and what stands to their right is a list of exponent arrays. The result of this 
action will be a list of True/False results with one result for each exponent array in 
the list. Below is a short Mathematica run that illustrates how the further modified 
function acts on lists. 

k = {4,5,l}; 

ja = {3,4,l}; jb = {1,2,3}; jc = {1,2,1}; 
list = {ja, jb, jcj; 

Apply[And,Thread[#l <= #2 & [#,k]]] & /@ list 

{True, False, True} (7.54) 

Observe that the output answer list is {True, False, True} because {3, 4, 1} < {4, 5, 1} 
is True, {1, 2, 3} < {4, 5, 1} is False, and {1, 2, 1} < {4, 5, 1} is True. 

• What we would really like to know is where the True items are in the list, because 
that will tell us where the j that satisfy j < k reside. This can be accomplished 
by use of the Mathematica command Position in conjunction with the result True. 
Below is a short Mathematica run that illustrates how this works. 

k = {4,5,l}; 

ja = {3,4,l};jb = {l,2,3};jc = {l,2,l}; 
list = {ja, jb, jc}; 

Position[Apply[And, Thread[#l <= #2 & [#,k]]] & /(§ list, True] 

{{1}, {3}} (7.55) 

Note that the output is an array of positions in the list for which j < k. There 
is, however, still one defect. Namely, the output array is an array of single-element 
subarrays, and we would like it to be simply an array of location numbers. This defect 
can be remedied by appending the Mathematica command Flatten, preceded by //, 
to the instruction string in (7.55). The short Mathematica run below illustrates this 
modification. 

k = {4,5,l}; 

ja= {3,4, 1}; jb = {1,2,3}; jc = {1,2,1}; 
list = {ja, jb, jc}; 

Position[Apply[And, Thread[#l <= #2 & [#, k]]] & /@ list, True]//Flatten 

{1, 3} (7.56) 

Now the output is a simple array containing the positions in the list for which j < k. 



59 



• The last step is to employ the ingredients in (7.56) to define the operator JSK[list, k]. 
The short Mathematica run below illustrates how this can be done. 



k = {4,5,1}; 

ja= {3,4, 1}; jb = {1,2,3}; jc = {1,2, l}; 
list = {ja, jb, jc}; 
JSK[list_, k_] := 

Position[Apply[And,Thread[#l <= #2 & [#,k]]] & /@ list, True]//Flatten; 
JSK[list, k] 



Lines 4 and 5 above define the operator JSK[list,k], line 6 invokes it, and line 7 
displays its output, which agrees with the output of (7.56). 

• With the operator JSK[list,k] in hand, we are prepared to construct tables B and 
Brev that will contain the B^ and the Brevk- The short Mathematica run below 
illustrates how this can be done. 



The first line employs the Mathematica command Table in combination with an 
implied Do loop to produce a two-dimensional array B. Values of k in the range [1, L] 
are selected sequentially. For each k value the associated exponent array k{k) = 
GAMMA[[k]] is obtained. The operator JSK then searches the full GAMMA array to find 
the list of r values associated with the j < k. All these r values are listed in a row. 
Thus, the array B consists of list of L rows, of varying width. The rows are labeled by 
k E [1, and in each row are the r values associated with the j < k. In the second 
line the Mathematica command Reverse is applied to B to produce a second array 
called Brev. Its rows are the reverse of those in B. For example, as the Mathem,atica 



{1, 3} 



(7.57) 



B = Table[JSK[GAMMA,GAMMA[[k]]], {k, 1,L, l}]; 

Brev = Reverse /@ B; 

B[[8]] 

Brev[[8]] 

B[[28]] 

Brev[[28]] 

{1,3,8} 

{8,3,1} 

{1,2,3,4,6,7,8,9,14,15,18,28} 
{28, 18,15,14,9,8,7,6,4,3,2,1} 



(7.58) 
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run illustrates, B[[8]], which is the 8*'* row of B, contains the list {1, 3, 8}, and Brev[[8]] 
contains the list {8, 3, 1}. Inspection of the r = 8 monomial in Table 4, that with 
exponents {0, 2, 0}, shows that it has the monomials with exponents {0,0,0}, {0,1,0}, 
and {0,2,0} as factors. And further inspection of Table 4 shows that the exponents 
of these factors have the r values {1,3,8}. Similarly B[[28]], which is the 28*'* row 
of B, contains the same entries that appear in the first column of Table 6. And 
Brev[[28]], which is the 28*'* row of Brev, contains the same entries that appear in 
the last column of Table 6. 

Finally, we need to explain how the arrays B and Brev can be employed to carry out 
polynomial multiplication. This can be done using the Mathematica dot product command: 

• The exhibit below shows a simple Mathematica run that illustrates the use of the dot 
product command. 

Unprotect[V]; 

U={.1,.2,.3,.4,.5,.6,.7,.8}; 

V = {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8}; 

U.V 

u = {l,3,5}; 
v = {6,4,2}; 
U[[n]] 
V[[v]] 

U[[u]].V[[v]] 
5.64 

{.1,.3,.5} 
{1.6,1.4,1.2} 

1.18 (7.59) 

As before, V must be unprotected. See line 1. The rest of the first part this run (lines 
2 through 4) defines two vectors U and V and then computes their dot product. Note 
that if we multiply the entries in U and V pairwise and add, we get the result 

.1 X 1.1 + .2 X 1.2 + • • • + .8 X 1.8 = 5.64, 

which agrees with the Mathematica result for U • V. See line 10. 

The second part of this Mathematica run, lines 5 through 9, illustrates a powerful 
feature of the Mathematica language. Suppose, as illustrated, we define two arrays u 
and V of integers, and use these arrays as arguments for the vectors by writing U[[u]] 
and V[[v]]. Then Mathematica uses the integers in the two arrays u and v as labels 
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to select the corresponding entries in U and V, and from these entries it makes new 
corresponding vectors. In this example, the 1^*, 3'''^, and S*'* entries in U are .1, .3, 
and .5. And the 6*'*, 4*^*, and 2"*^ entries in V are 1.6, 1.4, and 1.2. Consequently, we 
find that 

U[[u]] = {.1,.3,.5}, 

V[[v]] = {1.6,1.4,1.2}, 

in agreement with lines 11 and 12 of the Mathematica results. Correspondingly, we 
expect that U[[u]] • V[[v]] will have the value 

U[[u]] • V[[v]] = .1 X 1.6 + .3 X 1.4 + .5 X 1.2 = 1.18, 

in agreement with the last line of the Mathematica output. 

• Now suppose, as an example, that we set A; = 8 and use B[[k]] and Brev[[k]] in place 
of the arrays u and v. The Mathematica fragment below shows what happens when 
this is done. 



k = 8; 

B[[k]] 

Brev[[k]] 

U[[B[[k]]]] 

V[[Brev[[k]]]] 

U[[B[[k]]]].B[[Brev[[k]]]] 

{1,3,8} 

{8,3,1} 

{.1,.3,.8} 

{1.8,1.3,1.1} 

1.45 (7.60) 

From (7.58) we see that B[[8]] = {1,3,8} and Brev[[8]] = {8,3, 1} in agreement with 
lines 7 and 8 of the Mathematica output above. Also, the 1**, 3^*^, and 8*'^ entries 
in U are .1, .3, and .8. And the S*'', 3'"'^, and 1** entries in V are 1.8, 1.3, and 1.1. 
Therefore we expect the results 

U[[B[[k]]]] = {.l,.3,.8}, 

V[[Brev[[k]]]] ={1.8,1.3,1.1}, 
U[[B[[k]]]] • V[[Brev[[k]]]] = .1 X 1.8 + .3 x 1.3 + .8 x 1.1 = 1.45, 
in agreement with the last three lines of (7.60). 
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• Finally, suppose we carry out the operation U[[B[[k]]]] • V[[Brev[[k]]]] for all k G [l,L] 
and put the results together in a Tabic with entries labeled by k. According to 
(7.45), the result will be the pyramid for the product of the two polynomials whose 
individual pyramids are U and V. The Mathematica fragment below shows how this 
can be done to define a product function, called PROD, that acts on general pyramids 
U and V, using the command Table with an implied Do loop over A;. 

PROD[U_,V_] := Table [U[[B[[k]]]] • V[[Brev[[k]]]], {k, 1,L, 1}]; 
7.1.5 Implementation of Powers 

With operation of multiplication in hand, it is easy to implement the operation of raising a 
pyramid to a power. The code shown below in (7.61) demonstrates how this can be done. 

POWER[U_,0] := CI; 
POWER[U_, 1] := U; 
P0WER[U_,2] := PROD[U,U]; 
P0WER[U_,3] := PROD[U, POWER[U, 2]]; 

(7.61) 

Here CI is the pyramid for the Taylor series having one as its constant term and all other 
terms zero, 

01 = {1,0,0,0,- ••}. (7.62) 
It can be set up by the Mathematica code 

CI = Table [KroneckerDelta[k, l],{k, 1,L, 1}]; (7.63) 

which employs the Table command, the Kronecker delta function, and an implied Do loop 
over k. This code should be executed before executing (7.61), but after the value of L has 
been established. 



7.1.6 Replacement Rule 

Definition 2 The transformation A(z) A means replacement of every real variable Za 
in the arithmetic expression A[z) with an associated pyramid, and of every operation on 
real variables in A{z) with the associated operation on pyramids. 

Automatic Differentiation is based on the following corollary: if A{z) A, then A is the 
pyramid of A{z). 

For simplicity, we will begin our discussion of the replacement rule with examples 
involving only a single variable z. In this case monomial labeling, the relation between 
labels and exponents, is given directly by the simple rules 

r{j) = 1 + j and j(r) = r — 1. (7.64) 
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Table 7: A labeling scheme for monomials in one variable. 



r J 

1 

2 1 

3 2 

4 3 



See table 7. 

As a first example, consider the expression 

A = 2 + ^z* z). (7.65) 

We have agreed to consider the case m = 1. Suppose we also set p = 2, in which case 
L = 3. In ascending glex order, see Table 7, the pyramid for A is then 

2 + 3^^ A = (2,0,3). (7.66) 

Now imagine that A was not such a simple polynomial, but some complicated expression. 
Then the pyramid A could be generated by computing derivatives of ^ at z = and dividing 
them by the appropriate factorials. Automatic differentiation offers another way to find A. 
Assume that all operations in the arithmetic expression A have been encoded according 
to Definition 1. For our example, these are + and PROD. Let CI and Z be the pyramids 
associated with 1 and z, 

1 CI = (1,0,0), (7.67) 

2;-wZ = (0,1,0). (7.68) 

The quantity 2 + Sz^ results from performing various arithmetic operations on 1 and z. 
Definition 1 says that the pyramid of 2 + 3z^ is identical to the pyramid obtained by 
performing the same operations on the pyramids CI and Z. That is, suppose we replace 1 
and z with correctly associated pyramids CI and Z, and also replace * with PROD. Then, 
upon evaluating PROD, multiplying by the appropriate scalar coefficients, and summing, 
the result will be the same pyramid A, 

2 CI + 3 PROD[Z, Z] = A. (7.69) 

In this way, by knowing only the basic pyramids CI and Z (prepared beforehand), one can 
compute the pyramid of an arbitrary A{z). Finally, in contrast to numerical differentiation. 
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all numerical operations involved are accurate to machine precision. Mathematica code that 
implements (7.69) will be presented shortly in (7.70). 

Frequently, if A{z) is some complicated expression, the replacement rule will result 
in a long chain of nested pyramid operations. At every step in the chain, the pyramid 
resulting from the previous step will be combined with some other pyramid to produce a 
new pyramid. Each such operation has two arguments, and Definition 1 applies to each 
step in the chain. Upon evaluating all pyramid operations, the final result will be the 
pyramid of A{z). 

By using the replacement operator the above procedure can be represented as: 

1 CI, z ^ Z, A ^ k. 

The following general recipe then applies: In order to derive the pyramid associated with 
some arithmetic expression, apply the rule to all its variables, or parts, and replace all 
operations with operations on pyramids. Here "apply the rule" to something means 
replace that something with the associated pyramid. And the term "parts" means sub- 
expressions. Definition 1 guarantees that the result will be the same pyramid A no matter 
how we split the arithmetic expression A into sub-expressions. It is only necessary to 
recognize, in case of using sub-expressions, that one pyramid expression should be viewed 
function of another. 

For illustration, suppose we regard the A given by (7.65) to be the composition of two 
functions, F{z) = 2 + 32; and G{z) = z^, so that A{z) = F{G(z)). Instead of associating a 
constant and a single variable with their respective pyramids, let us now associate whole 
sub-expressions. In addition, let us label the pyramid expressions on the right of with 
with some names, F and G: 

2 + 3z-w2Cl + 3Z = F[Z] 
z'^ PROD[Z,Z] = G[Z] 

A{z) -w F[G[Z]] = A. 

We have indicated the explicit dependence on Z. It is important to note that F[Z] is 
a pyramid expression prior to executing any the pyramid operations, i.e it is not yet a 
pyramid, but is simply the result of formal replacements that follow the association rule. 
Mathematica code for the simple example (7.69) is shown below, 

CI = {1,0,0}; 

Z = {0,1,0}; 

2 CI + 3 PROD[Z,Z] 

{2,0,3} (7.70) 

Note that the result (7.70) agrees with (7.66). This example does not use any nested 
expressions. We will now illustrate how the same results can be obtained using nested 
expressions. 
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We begin by displaying a simple Mathematica program/execution, that employs ordi- 
nary variables, and uses Mathematica' s intrinsic abilities to handle nested expressions. The 
program/execution is 

f[z_] := 2 + 3z; 
g[z_] := z^; 

f[g[z]] 

2 + (7.71) 

With Mathematica the underscore in z_ indicates that z is a dummy variable name, and 
the symbols := indicate that f is defined with a delayed assignment. That is what is done 
in line one above. The same is done in line two for g. Line three requests evaluation of the 
nested function f{g{z)), and the result of this evaluation is displayed in line four. Note 
that the result agrees with (7.65). 

With this background, we are ready to examine a program with analogous nested 
pyramid operations. The same comments apply regarding the use of underscores and 
delayed assignments. The program is 

CI = {1,0,0}; 
Z = {0,1,0}; 
F[Z_] := 2 CI + 3 Z; 
G[Z_] := PROD[Z,Z]; 
F[G[Z]] 

{2,0,3} (7.72) 

Note that line (7.72) agrees with line (7.70), and is consistent with line (7.66). 

We close this subsection with an important consequence of the replacement rule and 
nested operations, which we call the Taylor rule. We begin by considering functions of a 
single variable. Suppose the function G{x) has the special form 

G{x) = z'^ + x (7.73) 

where z'^ is some constant. Let F be some other function. Consider the composite (nested) 
function A defined by 

A{x) = F{G{x)) = Fiz"^ + x). (7.74) 

Then, assuming the necessary analyticity, by the chain rule A evidently has a Taylor 
expansion in x about the origin of the form 

A = A{0) + A'{0)x + {l/2)A'\0)x'^ + --- 

= F(z'^) + F'{z'^)x + {l/2)F"{z'^)x^ + ■■■ . (7.75) 
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We conclude that if we know the Taylor expansion of A about the origin, then we also 
know the Taylor expansion of F about z*^, and vice versa. Suppose, for example, that 

F{z) = l + 2z + 3z'^ (7.76) 

and 

/ = 4. (7.77) 

Then there is the result 

A{x) = F{G{x)) = F{z'^ + x) = 1 + 2(4 + x) + 3(4 + xf = 57 + 26x + Sx"^. (7.78) 

We now show that this same result can be obtained using pyramids. The Mathematica 
fragment below illustrates how this can be done. 

CI = {1,0,0}; 
X = {0,1,0}; 
zd = 4; 

F[Z_] := 1 CI + 2 Z + 3 PROD[Z, Z]; 

G[X_] := zd CI + X; 

F[G[X]] 

{57,26,3} (7.79) 

Note that (7.79) agrees with (7.78). 

Let us also illustrate the Taylor rule in the two- variable case. Let F{zi,Z2) be some 
function of two variables. Introduce the functions G(xi) and H{xi) having the special 
forms 

G{xi) = zf + xi, (7.80) 
H{x2) = zi + X2, (7.81) 
where zf and are some constants. Consider the function A defined by 

A{xi,X2) = F{G{xi),H{x2)) = F{zf + xi, z| + X2). (7.82) 

Then, again assuming the necessary analyticity, by the chain rule A evidently has a Taylor 
expansion in xi and X2 about the origin (0, 0) of the form 

A = A{0,0) + [diA{0,0)]xi + [d2A{0,0)]x2 

+ (1/2)[(9i)2A(0, 0)]xl + [did2A{0, 0)]xiX2 + {l/2)[{d2fA{0, 0)]xl + • • • 
= F{zf, 4) + [OiF{zf, zi)]xi + [d2F{zf, zi)]x2 

+{l/2)[{difF{zf,z^)]xl + [did2AFizf,z^)]x,X2 + {l/2)[{d2fA{F{zf,zl))]xl + ■■■ 

(7.83) 
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where 



di = d/dxi, 82 = 8/ 8x2 



(7.84) 



when acting on A, and 



81 = 8/8zi, 82 = 8/8z2 



(7.85) 



when acting on F. We conclude that if we know the Taylor expansion of A about the origin 
(0, 0), then we also know the Taylor expansion of F about {zf, Zg), and vice versa. 
As a concrete example, suppose that 



Then, hand calculation shows that F{G{xi), H{x2)) takes the form 

F{zf + xi,zl + X2) = F{G{xi),H{x2)) 

= 899 + 98x1 +4x? + 134x2 + 5xia;2 + 6a;|. (7.88) 

Below is a Mathematica execution that finds the same result, 

F[zl_, z2_] := 1 + 2 zl + 3 z2 + 4 zl^ + 5 zl z2 + 6 z2^ 
G[xl_] := zdl +xl; 
H[x2_] := zd2 + x2; 
zdl = 7; 

zd2 = 8; 

A = F[G[xl],H[x2]] 
Expand [a] 

1 + 2 (7 + xl) + 4 (7 + xlf + 3 (8 + x2) + 5 (7 + xl) (8 + x2) + 6 (8 + x2)2 
899 + 98 xl + 4 xl^ + 134 x2 + 5 xl x2 + x2'^ 



F{zi,Z2) = 1 + 22;i + 32:2 + 4z? + 5ziZ2 + 6^1 



(7.86) 



and 



zf = 7, 4 = 8. 



(7.87) 



(7.89) 
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The calculation above dealt with the case of a function of two ordinary variables. We 
now illustrate, for the same example, that there is an analogous result for pyramids. For 
future reference. Table 8 shows our standard modified glex sequencing applied to the case 
of two variables. 



Table 8: A labeling scheme for monomials in two variables. 



r 


Ji 


J2 


1 








2 


1 





3 





1 


4 


2 





5 


1 


1 


6 





2 


7 


3 





8 


2 


1 


9 


1 


2 


10 





3 



Following the replacement rule, we should make the substitutions 

zf + xi-^zdl CI + XI, (7.90) 
zi + X2'-^ zd2 CI + X2, (7.91) 

1 + 2 zi + 3 Z2 + 4 + 5 zi Z2 + 6 -w 

CI + 2 Zl + 3 Z2 + 4 PR0D[Z1, Zl] + 5 PR0D[Z1, Z2] + 6 PR0D[Z2, Z2]. 

(7.92) 
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The Mathematica fragment below, executed for the case m = 2 and p = 2, in which case 
L = 6, illustrates how the analogous result is obtained using pyramids, 

CI = {1,0,0,0,0,0}; 
XI = {0, 1,0,0,0,0}; 
X2 = {0,0, 1,0,0,0}; 

F[Z1_, Z2_] := CI + 2 Zl + 3 Z2 + 4 PR0D[Z1, Zl] + 5 PR0D[Z1, Z2] 

+6 PR0D[Z2,Z2]; 

G[X1_] := zOl CI +X1; 

H[X2_] := z02 CI + X2; 

zdl = 7; 

zd2 = 8; 

F[G[X1],H[X2]] 

{899,98,134,4,5,6} (7.93) 

Note that, when use is made of Table 8, the last line of (7.93) agrees with (7.88) and the 
last line of (7.89). 

7.2 Replacement Rule and Numerical Integration 
7.2.1 Numerical Integration 

Consider the set of differential equations (1.5). A standard procedure for their numerical 
integration from an initial time V = t° to some final time is to divide the time axis into 
a large number of steps N, each of small duration h, thereby introducing successive times 
defined by the relation 

e = t° + nh with n = 0, 1, ■ • • , AT. (7.94) 

By construction, there will also be the relation 

Nh = tf - f. (7.95) 

The goal is to compute the vectors z", where 

z'^ = z(r), (7.96) 

starting from the vector 2;°. The vector z° is assumed given as a set of definite numbers, 
i.e. the initial conditions at t^. 

If we assume Poincare analyticity in t, we may convert the set of differential equations 
(1.5) into a set of recursion relations for the z" in such a way that the obtained by 
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solving the recursion relations differ from the true z" by only small truncation errors of 
order h"^. (Here m is not the number of variables, but rather some fixed integer describing 
the accuracy of the integration method.) One such procedure, fourth-order Runge Kutta 
(RK4), is the set of marching/recursion rules 



^n+i ^ ^ 1 ^ 26 + 2c + d) (7.97) 



where, at each step. 



a = hfiz^,n, (7.98) 
b = hfiz-+'^a,e + ^h), 

d = hf{z'^ + c,e + h). 

Thanks to the genius of Runge and Kutta, the relations (7.97) and (7.98) have been con- 
structed in such a way that the method is locally (at each step) correct through order /i^, 
and makes local truncation errors of order h^. 

In the case of a single variable, and therefore a single differential equation, the relations 
(7.97) and (7.98) may be encoded in the Mathematica form shown below. Here Zvar is the 
dependent variable, t is the time, Zt is a temporary variable, tt is a temporary time, and 
ns is the number of integration steps. The program employs a Do loop over i so that the 
operations (7.97) and (7.98) are carried out ns times. 
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RK4 := ( 
to = t; 

Do[ 



Aa 


= h F[Zvar, t]; 


Zt 


= Zvar + (l/2)Aa; 


tt 


= t + h/2; 


Bb 


= hF[Zt,tt]; 


Zt 


= Zvar + (l/2)Bb; 


Cc 


= hF[Zt,tt]; 


Zt 


= Zvar + Cc; 


tt 


= t + h; 


Dd 


= h F[Zt,ttl: 



Zvar = Zvar + (l/6)(Aa+ 2 Bb + 2 Cc + Dd); 

t = to + i h;, 

{i,l,ns,l} 

] 

) 

(7.99) 

7.2.2 Replacement Rule, Single Equation/ Variable Case 

We now make what, for our purposes, is a fundamental observation: The operations that 
occur in the Runge Kutta recursion rules (7.97) and (7.98) and realized in the code above 
can be extended to pyramids by application of the replacement rule. In particular, the 
dependent variable z can be replaced by a pyramid, and the various operations involved in 
the recursion rules can be replaced by pyramid operations. Indeed if wc look at the code 
above, apart from the evaluation of F, we see that the quantities Zvar, Zt, Aa, Bb, Cc, and 
Dd can be viewed, if we wish, as pyramids since the only operations involved are scalar 
multiplication and addition. The only requirement for a pyramidal interpretation of the 
RK4 Mathematica code is that the right side of the differential equation, F[*, *], be defined 
for pyramids. Finally, we remark that the features that make it possible to interpret the RK4 
Mathematica code either in terms of ordinary variables or pyramidal variables will hold for 
Mathematica realizations of many other familiar numerical integration methods including 
other forms of Runge Kutta, predictor-corrector methods, and extrapolation methods. 

To make these ideas concrete, and to understand their implications, let us begin with a 
simple example. Suppose, in the single variable case, that the right side of the differential 
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equation has the simple form 

f{z,t) = -2tz^. (7.100) 

The differential equation with this right side can be integrated analytically to yield the 
solution 

= z7[l + z°(t-t°)2]. (7.101) 
In particular, for the case i° = 0, 2;° = 1, and t = 1, there is the result 

z(l) = z7[l + z°] = l/2. (7.102) 

Let us also integrate the differential equation with the right side (7.100) numerically. 
Shown below is the result of running the associated Mathematica Runge Kutta code for 
this case. 



Clear["Global' * "]; 

F[Z_,t_] := -2 t Z^; 

h= .1; 

ns = 10; 

t = 0; 

Zvar = 1.; 

RK4; 

t 

Zvar 
1. 

0.500001 

(7.103) 

Note that the last line of (7.103) agrees with (7.102) save for a "1" in the last entry. As ex- 
pected, and as experimentation shows, this small difference, due to accumulated truncation 
error, becomes even smaller if h is decreased (and correspondingly, ns is increased). 

Suppose we expand the solution (7.102) about the design initial condition z*^^ = 1 by 
replacing z° by 2;'^° + x and expanding the result in a Taylor series in x about the point 
x=0. Below is a Mathematica run that performs this task. 

zdO = 1; 

Series[(zdO + x)/(l + zdO + x), {x, 0, 5}] 

1 X x^ x^ x^ 



2 + 4-y + 16 - 32 + 64 + ^[^l 



(7.104) 
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We will now see that the same Taylor series can be obtained by the operation of 
numerical integration applied to pyramids. The Mathematica code below shows, for our 
example differential equation, the application of numerical integration to pyramids. 

Clear["Global' * "]; 
Needs [" Combinatorica' "] ; 

m = 1; p = 5; 

GAMMA = Compositions [0,m]; 

Do[GAMMA = Join[GAMMA, Reverse[Compositions[d,m]]], {d, l,p, l}]; 
L = Length [GAMMA]; 
JSK[list_,k_] := 

Position[Apply[And,Thread[#l <= #2 & [#,k]]] & /@ list, True]//Flatten; 
B = Table [JSK[GAMMA, GAMMA[[r]]], {r, 1, L, 1}]; 
Brev = Reverse/® B; 

PROD[U_,V_] := Table [U[[B[[k]]]].V[[Brev[[k]]]],{k, 1,L, 1}]; 
F[Z_,t_] := -2 t PROD[Z,Z]; 
h= .01; 
ns = 100; 

t = 0; 
zdO = 1; 

CI = {1,0,0,0,0,0}; 
X = {0,1, 0,0, 0,0}; 
Zvar = zdO CI + X; 
RK4; 
t 

Zvar 
1. 

{0.5, 0.25, -0.125, 0.0625, -0.03125, 0.015625} (7.105) 

The first 11 lines of the code set up what should be by now the familiar procediue for 
labeling and multiplying pyramids. In particular, m = 1 because we are dealing with a 
single variable, and p = 5 since we wish to work through fifth order. The line 

F[Z_, t_] := -2 t PROD[Z, Z] (7.106) 

defines F[*, *] for the case of pyramids, and is the result of applying the replacement rule 
to the right side of / as given by (7.100), 

-2i ^ -2tPR0D[Z,Z]. (7.107) 
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Lines 13 through 15 are the same as hnes 3 through 5 in (7.103) except that, in order 
to improve numerical accuracy, the step size h has been decreased and correspondingly 
the number of steps ns has been increased. Lines 16 through 19 now initialize Zvar as a 
pyramid with constant part zdO and first-order monomial part 1, 

Zvar = zdO CI + X. (7.108) 

These lines are the pyramid equivalent of line 6 in (7.103). Finally lines 20 through 22 are 
the same as lines 7 through 9 in (7.103). In particular, the line RK4 in (7.103) and the line 
RK4 in (7.105) refer to exactly the same code, namely that in (7.99). 

Let us now compare the outputs of (7.103) and (7.105). Comparing the penultimate 
lines in each we see that the final time t = 1 is the same in each case. Comparing the 
last lines shows that the output Zvar for (7.105) is a pyramid whose first entry agrees 
with the last line of (7.103). Finally, all the entries in the pyramid output agree with the 
Taylor coefficients in the expansion (7.104). We see, in the case of numerical integration 
(of a single differential equation), that replacing the dependent variable by a pyramid, with 
the initial value of the pyramid given by (7.108), produces a Taylor expansion of the final 
condition in terms of the initial condition. 

What accounts for this near miraculous result? It's the Taylor rule described at the 
end of Section 7.1.6. We have already learned that to expand some function F{z) about 
some point z'^ we must evaluate F[z'^ + x). See (7.74). We know that the final Zvar, call 
it Zvar^"^, is an analytic function of the initial Zvar, call it Zvar^"^, so that we may write 

Zvar^"" = Zvar^^'iZvar'"') = g{Zvar''') (7.109) 

where g is the function that results from following the trajectory from t = to t = t^^. 

Therefore, by the Taylor rule, to expand Zvar^^ about Zvar™ = z'^'^, we must evaluate 
^^Q,j,fin^^do _|_ That, with the aid of pyramids, is what the code (7.105) accomplishes. 

7.2.3 Multi Equation/ Variable Case 

Because of Mathematica's built-in provisions for handling arrays, the work of the previous 
section can easily be extended to the case of several differential equations. Consider, as an 
example, the two- variable case for which / has the form 

fiiz,t) = -zl 

f2{z,t) = +2ziZ2. (7.110) 

The differential equations associated with this / can be solved in closed form to yield, with 
the understanding that t° = 0, the solution 

zi{t) = zl/{l + tzl), 

Z2{t) = zl{l + tz\f. (7.111) 
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For the final time t = 



1 we find the result 



zi(l) = z?/(l + ^?), 

Z2(1) = Z2°(1 + ^?)'- (7.112) 
Let us expand the solution (7.112) about the design initial conditions 

zf = 1, 

zf = 2, (7.113) 



by writing 



z^ = zf + xi = 1 + xi, 

z^ = zf + X2 = 2 + X2. (7.114) 

Doing so gives the results 

zi(l) = (l + xi)/(2 + xi) = (2 + xi-l)/(2 + xi) = l-l/(2 + xi) = 

= 1 - (1/2)(1 + xi/2)-i = 1 - (1/2)[1 - xi/2 + (xi/2)2 - (xi/2)3 + • • • 
= (1/2) + il/4)xi - (l/8)x? + (l/16)x? + • • • , 

(7.115) 



Z2{1) = (2 + a;2)(2 + xi)2 

= 8 + 8xi+4x2 + 2xi + 4xiar2 + a;?a;2. (7.116) 

We will now explore how this same result can be obtained using the replacement rule 
applied to the operation of numerical integration. As before, we will label individual 
monomials by an integer r. Recall that Table 8 shows our standard modified glex sequencing 
applied to the case of two variables. 

The Mathematica code below shows, for our two- variable example differential equation, 
the application of numerical integration to pyramids. Before describing the code in some 
detail, we take note of the bottom two lines. When interpreted with the aid of Table 8, 
we see that the penultimate line of (7.117) agrees with (7.115), and the last line of (7.117) 
nearly agrees with (7.116). The only discrepancy is that for the monomial with label r = 7 
in the last line of (7.117). In the Mathematica output it has the value —1.16563x10"^ while, 
according to (7.116), the true value should be zero. This small discrepancy arises from the 
truncation error inherent in the RK4 algorithm, and becomes smaller as the step size h 
is decreased (and ns is correspondingly increased), or if some more accurate integration 
algorithm is used. We conclude that, with the use of pyramids, it is also possible in the 
two-variable case to obtain Taylor expansions of the final conditions in terms of the initial 
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conditions. Indeed, what is involved is again the Taylor rule applied, in this instance, to 
the case of two variables. 

Clear["Global' * "]; 
Needs["Combinatorica' "]; 
m = 2; p = 3; 

GAMMA = Compositions[0,m]; 

Do[GAMMA = Join[GAMMA,Reverse[Compositions[d,m]]], {d, l,p, 1}]; 
L = Length [GAMMA]; 

JSK[list_,k_] : = 

Position[Apply[And, Thread[#l <= #2 & [#, k]]] & /@ list, True]//Flatten; 
B = Table[JSK[GAMMA,GAMMA[[r]]],{r, 1,L, l}]; 
Brev = Reverse/® B; 

PROD[U_,V_] :=Table[U[[B[[k]]]].V[[Brev[[k]]]],{k,l,L,l}]; 
F[Z_,t_] := {-PR0D[Z[[1]],Z[[1]]],2. PROD[Z[[l]], Z[[2]]]}; 
h = .01; 
ns = 100; 
t = 0; 

zdO = {l.,2.}; 

CI = Table [KroneckerDelta[k, l], {k, 1, L, l}]; 
X[l] = Table[KroneckerDelta[k, 2], {k, 1, L, l}]; 
X[2] = Table[KroneckerDelta[k, 3], {k, 1, L, l}]; 
Zvar = {zdO[[l]] CI + X[l], zd0[[2]] CI + X[2]}; 
RK4; 
t 

Zvar 
1. 

{{0.5, 0.25, 0., -0.125, 0., 0., 0.0625, 0., 0., 0, }, 
{8., 8., 4., 2., 4., 0., -1.16563 x 10"^ 1., 0., 0.}} (7.117) 

Let us compare the structures of the routines for the single variable case and multi 
(two) variable case as illustrated in (7.105) and (7.117). The first difference occurs at line 
3 where the number of variables m and the maximum degree p are specified. In (7.117) 
m is set to 2 because we wish to treat the case of two variables, and p is set to 3 simply 
to limit the lengths of the output arrays. The next difference occurs in line 12 where the 
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right side F of the differential equation is specified. The major feature of the definition of 
F in (7.117) is that it is specified as two pyramids because the right side of the definition 
has the structure {*,*} where each item * is an instruction for computing a pyramid. In 
particular, the two pyramids are those for the two components of / as given by (7.110) 
and use of the replacement rule, 



The next differences occur in lines 16 through 20 of (7.117). In line 16, since specification 
of the initial conditions now requires two numbers, sec (7.113), zdO is specified as a two- 
component array. In lines 17 and 18 of (7.105) the pyramids CI and X are set up explicitly 
for the case p = 5. By contrast, in lines 17 through 19 of (7.117), the pyramids CI, X[l], 
and X[2] are set up for general p with the aid of the Table command and the Kroneckcr 
delta function. Recall (7.62) and observe from Tables 4, 7, and 8 that, no matter what the 
values of m and p, the constant monomial has the label r = 1 and the monomial xi has 
the label r = 2. Moreover, as long as m > 2 and no matter what the value of p, the X2 
monomial has the label r = 3. Finally, compare line 19 in (7.105) with line 20 in (7.117), 
both of which define the initial Zvar. We see that the difference is that in (7.105) Zvar is 
defined as a single pyramid while in (7.117) it is defined as a pair of pyramids of the form 
{*, *}. Most remarkably, all other corresponding lines in (7.105) and (7.117) are the same. 
In particular, the same RK4 code, namely that given by (7.99), is used in the scalar case 
(7.103), the single pyramid case (7.105), and the two-pyramid case (7.117). This multi-use 
is possible because of the convenient way in which Mathematica handles arrays. 

We conclude that the pattern for the multivariable case is now clear. Only the following 
items need to be specified in an m dependent way: 

• The value of m. 

• The entries in F with entries entered as an array {*, *,•••} of m pyramids. 

• The design initial condition array zdO. 

• The pyramids for CI and X[l] through X[m]. 

• The entries for the initial Zvar specified as an array 

{zdO[[l]] Cl-FX[l],zdO[[2]] Cl + X[2],--- ,zdO[[m]] Cl-FX[m]} of m pyramids. 

7.3 Duffing Equation Application 

Let us now apply the methods just developed to the case of the Duffing equation with 
parameter dependence as described by the relations (6.12) through (6.17). Mathematica 
code for this purpose is shown below. By looking at the final lines that result from executing 



-zf ^-PR0D[Z[[1]],Z[[1]]], 



(7.118) 



2ziZ2 2. PR0D[Z[[1]],Z[[2]]]. 



(7.119) 
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this code, we see that the final output is an array of the form {{*}, {*}, {*}}• That is, the 
final output is an array of three pyramids. This is what we expect, because now we are 
dealing witli three variables. See line 3 of the code, which sets m=3. Also, for convenience 
of viewing, results are calculated and displayed only through third order as a consequence 
of setting p = 3. 

Clear["Global' * "]; 
Needs["Combinatorica' "]; 
m = 3; p = 3; 

GAMMA = Compositions[0,m]; 

Do[GAMMA = Join[GAMMA,Reverse[Compositions[d,m]]], {d, l,p, 1}]; 
L = Length [GAMMA]; 
JSK[list_,k_] := 

Position[Apply[And, Thread[#l <= #2 & [#, k]]] & /@ list, True]//Flatten; 
B = Table [JSK [GAMMA, GAMMA[[r]]], {r, 1, L, 1}]; 
Brev = Reverse/® B; 

PROD[U_,V_] := Table[U[[B[[k]]]].V[[Brev[[k]]]],{k, 1,L, l}]; 
POWER[U_, 2] := PROD[U,U]; 
POWER[U_, 3] := PROD[U, POWER[U, 2]]; 
CO = Table [0,{k, 1,L, l}]; 
F[Z_,t_] := {Z[[2]], 

-2. betaPR0D[Z[[3]],Z[[2]]] - PR0D[P0WER[Z[[3]], 2], Z[[l]]]- 

P0WER[Z[[1]],3] - eps Sin[t] P0WER[Z[[3]], 3], 

CO}; 
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ns = 100; 
t = 0; 

h = (2Pi)/ns; 

beta = .1; eps = 1.5; 

zdO = {.3, .4, .5}; 

CI = Table[KroneckerDelta[k, l],{k, 1,L, l}]; 

X[l] = Table[KroiieckerDelta[k, 2],{k, 1,L, 1}]; 

X[2] = Table[KroneckerDelta[k,3],{k, 1,L, l}]; 

X[3] = Table[KroneckerDelta[k,4],{k, 1,L, l}]; 

Zvar = {zdO[[l]] CI + X[l], zd0[[2]] CI + X[2], zd0[[3]] CI + X[3]}; 

RK4; 

t 

Zvar 

27r 

{{-0.0493158, 0.973942, -0.110494, 5.51271, 3.54684, 3.46678, 
11.2762, 2.36463, 1.0985, 23.3332, -1.03541, -3.23761, -12.8064, 
4.03421, -23.4342, -17.8967, 1.96148, 5.07403, -36.9009, 25.1379}, 
{0.439713, 1.05904, 0.427613, 3.3177, 0.0872459, 0.635397, -3.02822, 

1.77416, -4.10115, 3.16981, -2.43002, -5.33643, -7.77038, -6.08476, 
- 0.541465, -21.1672, -1.4091, -9.54326, 14.6334, -39.2312}, 
{0.5, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}} 

(7.120) 

The first unusual fragments in tlie code are lines 12 and 13, which define functions that 
implement the calculation of second and third powers of pyramids. Recall Section 7.1.5. 
The first new fragment is line 14, which defines the pyramid CO with the aid of the Table 
command and an implied Do loop. As a result of executing this code, CO is an array of 
L zeroes. The next three lines, lines 15 through 18, define F, which specifies the right 
sides of equations (6.12) through (6.14). Sec (6.15) through (6.18). The right side of F is 
of the form {*,*,*}, an array of three pyramids. By looking at (6.15) and recalling the 
replacement rule, we see that the first pyramid should be Z[[2]], 

^2^Z[[2]]. (7.121) 

The second pyramid on the right side of F is more complicated. It arises by applying the 
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replacement rule to the right side of (6.16) to obtain the associated pyramid, 



— 2/3z3Z2 — z'^zi — zf — ez| sint 

-2. betaPR0D[Z[[3]],Z[[2]]] - PR0D[P0WER[Z[[3]], 2], Z[[l]]]- 
P0WER[Z[[1]], 3] - eps Sin[t] P0WER[Z[[3]], 3]. (7.122) 

The third pyramid on the right side of F is simplicity itself. From (6.17) wc sec that this 
pyramid should be the result of applying the replacement rule to the number 0. Hence, 
this pyramid is CO, 

0-^C0 = {0,0,--- ,0}. (7.123) 

The remaining lines of the code require little comment. Line 20 sets the initial time to 

0, and line 21 defines h in such a way that the final value of t will be 2tt. Line 22 establishes 
the parameter values f3 = .1 and e = 1.5, which are those for Figure 4. Line 23 specifies 
that the design initial condition is 

zi{0) = zf = .3, Z2{0) = zf = .4, Z3(0) = zf = .5 = a, (7.124) 

and consequently 

oj = l/a = 2. (7.125) 
See (6.3). Also, it follows from (6.2) and (6.5) that 

^(0) = ujQiO) = ojzi{0) = (2)(.3) = .6, (7.126) 

g'(0) = u;2g(0) = a;2^2(0) = (2^)(.4) = 1.6. (7.127) 

Next, lines 24 through 28 specify that the expansion is to be carried out about the initial 
conditions (7.124). Finally, line 29 invokes the RK4 code given by (7.99). That is, as before, 
no modifications are required in the integration code. 

A few more comments about the output are appropriate. Line 32 shows that the final 
time t is indeed 27r, as desired. The remaining output lines display the three pyramids that 
specify the final value of Zvar. From the first entry in each pyramid we see that 

zi(27r) = -0.0493158, (7.128) 

Z2{2tt) = 0.439713, (7.129) 

Z3(27r) = .5, (7.130) 

when there are no deviations in the initial conditions. The remaining entries in the pyra- 
mids are the coefficients in the Taylor series that describe the changes in the final conditions 
that occur when changes are made in the initial conditions (including the parameter a). 

We are, of course, particularly interested in the first two pyramids. The third pyramid has 
entries only in the first place and the fourth place, and these entries are the same as those 
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in the third pyramid pyramid for Zvar at the start of the integration, namely those in 
zd0[3] CI + X[3]. The fact that the third pyramid in Zvar remains constant is the expected 
consequence of (6.17). 

We should also describe how the Ais employed in Section 6.2 was actually computed. 
It could have been computed by setting p = 8 m (7.120) and specifying a large number of 
steps ns to insure good accuracy. Of course, when p = 8, the pyramids are large. Therefore, 
one does not usually print them out, but rather writes them to files or sends them directly 
to other programs for further use. 

However, rather than using RK4 in (7.120), we replaced it with an adaptive 4-5*^ order 
Runge-Kutta-Fehlberg routine that dynamically adjusts the time step h during the course 
of integration to achieve a specified local accuracy, and we required that the error at each 
step be no larger than 10~^^. Like the RK4 routine, the Runge-Kutta-Fehlberg routine, 
when implemented in Mathematica, has the property that it can integrate any number of 
equations both in scalar variable and pyramid form without any changes in the code.^'^ 

7.4 Relation to the Complete Variational Equations 

At this point it may not be obvious to the reader that the use of pyramids in integration 
routines to obtain Taylor expansions is the same as integrating the complete variational 
equations. We now show that the integration of pyramid equations is equivalent to the 
forward integration of the complete variational equations. For simplicity, we will examine 
the single variable case with no parameter dependence. The reader who has mastered this 
case should be able to generalize the results obtained to the general case. 
In the single variable case with no parameter dependence (2.1) becomes 



Let z'^{t) be some design solution and introduce a deviation variable C by writing 



z = f{z,t). 




z = z'^ + C- 



(7. 



132) 




(7. 



133) 



Also, the relations (2.4) and (2.5) take the form 




(7. 



134) 



where g has an expansion of the form 



oo 




(7. 



135) 



A Mathematica version of this code is available from the first author upon request. 
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Finally, (2.6) and (2.7) become 

z'' = f{z'',t), (7.136) 

oo 

C = giz'',t,0 = Y,9'it)e, (7.137) 

and (2.8) becomes 

oo 

c = Y.h'{tmy. (7.138) 

Insertion of (7.138) into both sides of (7.137) and equating like powers of Q now yields the 
set of differential equations 

oo 

h'"{t) = Y,9'mfm withi,/ > 1 (7.139) 
where the (universal) functions C/^ {h^) are given by the relations 



yji=i ) j"=i 

The equations (7.136) and (7.139) are to be integrated from t = t^^ = to t = i^" with 
the initial conditions 

z'^it^) = z'^, (7.141) 

/i\i°) = 1, (7.142) 

/i^"(t°) = Ofor/ > 1. (7.143) 

Let us now consider the numerical integration of pyramids. Upon some reflection, we see 
that the numerical integration of pyramids is equivalent to finding the numerical solution 
to a differential equation with pyramid arguments. For example, in the single- variable case, 
let Zvar(t) be the pyramid appearing in the integration process. Then, its integration is 
equivalent to solving numerically the pyramid differential equation 

{d/dt)Zvax{t) = F(Zvar, t). (7.144) 

We now work out the consequences of this observation. By the inverse of the replace- 
ment rule, we may associate a Taylor series with the pyramid Zvar(t) by writing 

Zvar(i) co{t) + E Cjit)x^. (7.145) 
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By (1.45) it is intended that the entries in the pyramid Zvar(i) be used to construct 
a corresponding Taylor series with variable x. In view of (7.108), there are the initial 
conditions 

co(to) = z'^ito), (7.146) 

ci(to) = 1, (7.147) 

Cj(to) = for J > 1. (7.148) 

Wc next seek the differential equations that determine the time evolution of the Cj{t). 
Under the inverse replacement rule, there is also the correspondence 

(d/dt)Zvar(i) co(t) + Cj{t)x^. (7.149) 

We have found a representation for the left side of (7.144). We need to do the same for the 
right side. That is, we need the Taylor series associated with the pyramid F(Zvar, i). By 
the inverse replacement rule, it will be given by the relation 

F(Zvar,t) ^ f{Y,Cj{t)x^,t). (7.150) 

j>0 

Here it is understood that the right side of (7.150) is to be expanded in a Taylor series 
about X = 0. From (7.134), (7.135), and (7.140) we have the relations 

f{Ycj{t)x^,t) = f{co{t))+g{co{t),t,Y,Cj{t)x^) 

j>0 j>i 

= fico{t)) + Y,9Ht)Q2cj{t)x^))'' 
k>i j>i 

k>i j>i 



(7.151) 



Therefore, there is the inverse replacement rule 



F(Zvar,i) - /(co(i)) + ^ /(i) ^ C/^(q)x^' (7.152) 

k>l j>l 

Upon comparing like powers of x in (7.149) and (7.152), we see that the pyramid 
differential equation (7.144) is equivalent to the set of differential equations 

co{t) = /(co(0), (7.153) 

Cj{t) = J29Hmi{ce). (7.154) 

ifc>i 
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Finally, compare the initial conditions (7.141) through (7.143) with the initial conditions 
(7.146) through (7.148), and compare the differential equations (7.136) and (7.139) with 
the differential equations (7.153) and (7.154). We conclude that that there must be the 
relations 



We have verified, in the single variable case, that the use of pyramids in integration routines 
is equivalent to the solution of the complete variational equations using forward integration. 
As stated earlier, verification of the analogous m-variable result is left to the reader. We 
also observe the wonderful convenience that, when pyramid operations arc implemented 
and employed, it is not necessary to explicitly work out the forcing terms ga{t) and the 
functions {h^), nor is it necessary to set up the equations (3.6). All these complications 
are handled implicitly and automatically by the pyramid routines. 

8 Concluding Summary 

Poincare analyticity implies that transfer maps arising from ordinary differential equations 
can be expanded as Taylor series in the initial conditions and also in whatever parameters 
may be present. Section 2 showed that the determination of these expansions is equiva- 
lent to solving the complete variational equations, and Sections 3 and 4 showed that the 
complete variational equations can be solved either by forward or backward integration. 
Sections 5 and 6 applied this procedure for the Duffing stroboscopic map and found, re- 
markably, that an 8**^ order polynomial approximation to this map produced an infinite 
period doubling cascade and apparent strange attractor that closely resembled those of the 
exact map. A final section described computer methods for automatically setting up and 
numerically integrating the complete variational equations. 
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